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FOREWORD 


The SPACE-A program described in this report is a modified version 
of the trajectory and observation generation portion of the Sequential Position 
and Covariance Estimation (SPACE) program originally acquired from NASA. 
This document contains a detailed analysis of the equations and models used 
by the program, a brief description of the routines used in the program, as 
well as instructions for its operation. The principal applications of SPACE-A 
are the prediction of space-vehicle trajectories and the generation of observa¬ 
tional data plus other trajectory related information. 

There have been numerous changes made in the original program in 
addition to its adaptation to the IBM 7030 and its subsequent debugging. The 
modified SPACE-A program is the result of the efforts of the authors and 
L. E. Wilkie. S. Schwartz aided in the debugging and rewriting of certain 
routines. Special thanks and appreciation are due to R. K. Squires, D. S. 
Woolston, and other members of the Special Projects Branch, Theoretical 
Division of the Goddard Space Flight Center from whom the SPACE program 
was obtained. 

The work reported in this document was performed by The MITRE 
Corporation, Bedford, Massachusetts, for the Directorate of Planning and 
Technology, Electronic Systems Division, of the Air Force Systems Command 
under Contract AF 19(628)-5165. 


REVIEW AND APPROVAL 


Publication of this technical report does not constitute Air Force approval of 
the report's findings or conclusions. It is published only for the exchange and 
stimulation of ideas. 

ANTHONY P. TRUNFIO 

Technical Advisor, Development Engineering Div. 

Directorate of Planning and Technology 
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ABSTRACT 

This report describes the SPACE-A program presently available for 
operational use. SPACE-A is the trajectory and observation generation 
portion of the Sequential Position and Covariance Estimation Program (SPACE) 
which is currently under development. The document contains a detailed 
analysis of the equations and models used by the program, a brief description 
of routines used in the program, as well as instructions for its operations. 

The principal applications of SPACE-A are the prediction of the space- 
vehicle trajectories and the generation of observational data plus other 
trajectory related information. 
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SECTION I 
OBJECTIVES 


INTRODUCTION 

The SPACE trajectory determination program was originally devel¬ 
oped for the Special Projects Branch* Theoretical Division, Goddard 
Space Flight Center by the Sperry Rand Systems Group, It was designed 
as a comprehensive trajectory determination and tracking program for 
both orbital and deep space flights. The original program consisted 
of three major modes of operation: 

(1) SPACE-A, designed for trajectory and observation generation 
without statistical processing, 

(2) SPACE-B1, designed for statistical estimation of the six 
state variables describing the position and velocity at 
prescribed points of the trajectory, 

(3) SPACE-B2, designed for statistical estimation of the six 
state variables describing position and velocity plus up 
to 20 additional states selected from a group of vari¬ 
ables which permit estimation of dynamic biases (para¬ 
meters affecting orbital motion) or observational biases 
(affecting measurements). 

The version of the SPACE program received by MITRE contained 
numerous errors in programming, SPACE-B1 and SPACE-B2 had never been 
debugged and many of the routines of SPACE-A were unworkable. In ad¬ 
dition, the programming had to be made compatible with the IBM 7030, 
Therefore, a large effort was put into checking the theory, debugging 
the program, and revising or rewriting certain routines. Many check¬ 
out runs were necessary during and after the debugging of the program. 

Since the SPACE-B2 mode includes the capabilities of the SPACE-B1 
mode, the SPACE-B1 mode was not debugged. The original SPACE-B2 
(which was not operational when received by MITRE) employed two dif¬ 
ferent statistical estimation techniques: 

(1) a minimum variance sequential filtering technique, the so- 
called "Kalman filter", and 
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(2) a non-recursive batch processing technique. 

Because a number of trajectory estimation programs at MITRE 
already employ the batch processing technique, the second option of the 
SPACE-B2 mode was left in a non-operational status. The first option 
of SPACE-B2, which incorporates the Kalman filter in the trajectory 
estimation procedure, is presently being debugged and checked out. 

The theory, programming, and results of its operation wil] be pub¬ 
lished in a succeeding document. The theoretical development and 
much of the programming of the trajectory generation portion of 
SPACE-B2 is identical with that used in SPACE-A. 

This document deals with the modified SPACE-A trajectory genera- 
tion program written in Fortran IV and is presently operational on 
the IBM 7030 at MITRE. This report is based on the original documents^ 5 
published by Sperry Rand for NASA. 

Section I contains a brief discussion of the objectives and capa¬ 
bilities of SPACE-A. Section II discusses the theoretical background, 
as well as the equations and methods employed by the program Section 
III consists of a user’s guide for operating the program. Section TV 
gives the program structure and a brief description of the function 
of each subroutine. 

The larger part of the report is contained in Section II since 
understanding or following the programming is really a matter of fol¬ 
lowing the appropriate equations. Although parts of Section II arc 
original, the contents of pertinent sections of the SPACE Analytic 
Manual L ^ have been freely used or modified to fit the description 
of the present program. Sections III and IV closely follow the format 
of the original User’s Manual and Programmer’s Manual however, 

since a number of changes in programming have been made, only Section 
III and Section IV of this report should be used in operating the modi¬ 
fied SPACE-A trajectory generation program. 
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CAPABILITIES 


The modified SPACE-A trajectory generation may be run in any of 
the following modes: 

(1) Trajectory generation - the computation of position and 
velocity of the vehicle at regular time points along a 
trajectory, 

(2) Observation generation - the computation of certain 
ground-based observations of the vehicle at regular 
time points, 

(3) Simulated data - the computation and writing onto 
tape of the observations in a format suitable for 
processing data by SPACE-B2. 

(4) Visibility computation - the time of initial view 
of the vehicle by a given station, the total time 
in sight by the station, and the observations 
while the vehicle is in sight. 

The required input consists of the initial date and initial time 
of a particular run and the initial position and velocity of the vehi¬ 
cle, When the effect of drag is to be included which is usually the 
case, the effective area and mass of the vehicle and the magnitude of 
solar flux must be specified. When the effect of direct solar radia¬ 
tion is included, an effective area pertaining to solar radiation must 
be specified. For observation calculations, the specification of the 
station location on the surface of the Earth is necessary. The de¬ 
tails of the units and format of these required inputs along with 
other optional input features are discussed in Section III, The types 
of output data, which consist primarily of trajectory and/or observa¬ 
tion information, is also given in Section III, 

The philosophy behind the modified SPACE-A trajectory generation 
program is that its primary use would be for Earth orbital trajec¬ 
tories, Toward this objective, its dynamic model includes the effects 
of the primary terra and oblateness perturbations of the Earth's grav¬ 
ity, the solar, lunar, and planetary gravitational perturbations, the 
effect of direct solar radiation and drag, as well as such dynamic ef¬ 
fects as precession, nutation, and the daily rotation of the Earth, 
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The model used for drag is a combination of the U.S. Standard 
Atmosphere 1962 and the Harris-Priester model for the upper atmosphere, 
which depends on the magnitude of solar flux* The model for Earth 
gravitational oblateness may include up to nine zonal harmonics, four 
sectorial harmonics, and twelve additional tesseral harmonics* The 
types of ground station observations that may be specified are given 
in Section II (see OBSERVATIONS). 

Although the modified SPACE-A program is intended primarily for 
Earth orbital trajectories, it may be used for other types of missions. 
However, a number of options that were included in the original SPACE-A 
trajectory generation program are not presently operational. These op¬ 
tions were a capability for modeling the perturbations due to thrust, 
a model for lunar libration, simple models for lunar gravitational ob¬ 
lateness perturbations, as well as models for the drag of Mars and 
Venus, and a capability for computing on-board observations from a 
vehicle. Although these options are not presently available in the 
modified SPACE-A program, with some effort they could also be included. 
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SECTION II 
THEORY 


GENERAL DESCRIPTION OF TRAJECTORY COMPUTATION 

Orbit prediction or trajectory computation is the process of cal¬ 
culating the position and velocity of a vehicle at any time subsequent 
to some initial time, provided the initial position and velocity of 
the spacecraft are given. 

To accomplish this prediction, one uses the laws of celestial 
mechanics embodied in the differential equations of motion. The 
forcing functions for the equations of motion are obtained from a dy¬ 
namic model which accounts for the accelerations on the vehicle. A 
reference frame is erected to express the components of the various 
vector quantities, and the equations of motion are numerically inte¬ 
grated subject to the given initial conditions. Once the vehicle tra¬ 
jectory is determined, the program can then generate observations. 

The coordinate system used in this program is based upon the 
position of the mean Earth’s equator and equinox at O^ 1 1 January of 

the year subsequent to the initial input time of the program. Coordi¬ 
nate directions of this frame are inertial with respect to the fixed 
stars; the center of origin of the system, however, may be transferred 
from one central body to another so that the spacecraft motion is 
specified relative to a point mass which itself has a proper motion. 
This reference frame is called the Base Date System or simply the in¬ 
ertial system. 

Observations made from the Earth are necessarily in a system 
different from the Base Date System. To accomplish transfer between 
various coordinate systems, transformations are provided (see COORDI¬ 
NATE SYSTEMS AND TRANSFORMATIONS). These transformations include the 
dynamical effects of precession, nutation, and daily rotation of the 
Earth. 
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All accelerations acting on the vehicle are specified in the Base 
Date System. The gravitational attractions of bodies in the solar sys¬ 
tem are functions only of position with respect to the vehicle; conse¬ 
quently, the program employs an ephemeris giving planetary coordinates 
relative to the Sun, and lunar coordinates relative to the Earth, all 
in a Base Date System. A Base Date System is specified for overlapping 
two-year blocks of data, the date corresponding to the middle of the 
two-year file. Specifying an initial time causes the program to choose 
an ephemeris file having as its Base Date the beginning of the year 
following the initial input time of the program. In this way, at least 
one full year of ephemeris information is available before a change of 
reference system is necessary. 

Another acceleration specified in the Base Date System without 
transformation is that arising from solar radiation pressure. Since 
this acceleration is a function of relative position between the Sun 
and the vehicle, its direction is given in the proper frame by using 
information from the ephemeris. 

Other accelerations, such as perturbations due to Earth oblate¬ 
ness effects must be transformed through nutation and precession to the 
proper frame. All positions, velocities and accelerations are computed 
in canonical units (i.e.. Earth radii (ER), Earth radii per hour 
(ER/hr), Earth radii per hour per hour (ER/hr 2 ), respectively) and ap¬ 
propriate constants are used to transform to other units of measure. 

Vehicle motion is always computed relative to some reference body; 
a planet; the Moon; the Sun. Consequently, the equations of motion 
contain a term which accounts for the acceleration of the reference 
body on the vehicle. The remaining accelerations are usually, but not 
always, much smaller than this primary acceleration and are therefore 
called perturbations. In most cases, they can be regarded as giving 
rise to small disturbances in the orbit determined by the reference 
body and the initial conditions. 
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Reference bodies are changed during a trajectory calculation when 
the vehicle leaves the "region of influence" associated with a partic¬ 
ular body. Regions of influence are computed for a body with respect 
to the object of which it is a satellite. Hence, each planet has a 
region of influence defined relative to the Sun, and the Moon h^s a 
similar region defined relative to the Earth, In transferring into or 
out of such a region, velocity as well as position with respect to the 
new reference body must be calculated. 

Since no analytic solution exists for the equations of motion, 
numerical methods are employed to compute the components of position 
and velocity. In the program, a choice may be made between using 
straightforward integration and using Encke's method. The former 
technique, called Cowell's method, is conceptually simple, but suffers 
from precision and machine running time problems. Encke’s method, 
although somewhat more complicated, gives dividends in both precision 
and machine efficiency. In this procedure, the Keplerian orbit 
arising trom the reference body central force field is taken as a 
nominal trajectory. The perturbation accelerations are integrated 
and the resulting position and velocity increments are added to the 
Keplerian solutions. Naturally, Encke's method is most effective when 
the perturbations are small. 

The present program has the capability to compute a number of 
ground-based observation::. Corrections are provided for the refraction 
of an electromagnetic signal by the troposphere and by the ionosphere. 
Adjustments are computed for errors in elevation angle, range, and 
radial velocity; other angular corrections are calculated from the 
adjustment in elevation angle. 

Description of the basic equations used by the program for the 
dynamic model now follows. It includes derivations of the accelera¬ 
tions due to the planets, Earth oblateness, drag, and direct solar 
radiation pressure. It also discusses the coordinate systems, numeri¬ 
cal integration methods, and observation calculations including the 
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refraction model used by the program* 

DYNAMIC MODEL 

The equations of motion for a space vehicle are second-order 
differential equations which describe the accelerations arising from 
the forces acting on the vehicle. These forces can be classified as 
follows: 

(1) Gravitational acceleration of the reference body - primary. 

(2) Gravitational perturbations due to other bodies (e.g., 
planets). 

(3) Gravitational perturbations due to reference body 
oblateness. 

(4) Perturbations due to drag. 

(5) Perturbations due to direct solar radiation pressure 
on the space vehicle. 

For orbit determination the reference body primary gravitational 
force almost always predominates the perturbation forces given above. 
The only exception occurs during reentry of a vehicle into the atmos¬ 
phere, when drag forces may be larger than the primary force of gra¬ 
vity. Another force which may exceed the primary force of gravity is 
that of vehicle thrust. Since this program is designed primarily for 
orbital computations, this force is not discussed here although it 
can be readily included if desired. 

The simplest gravitational force field is that due to a single 
point mass or equivalently that due to a homogeneous ponderable body 
which is perfectly spherical. In this case the equations of motion 
of a vehicle with respect to the ponderable reference body are: 


where 

y 

G 

M 

m 



G(M+m) = GM, since m << M 
is the universal gravitational constant, 
is the mass of the reference body, 
is the mass of vehicle. 
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R is the vector position of the vehicle with respect to 
the reference body center. 

With initial conditions Rq and Rq Equation (1) defines a "two- 
body" or Keplerian orbit, which arises from the primary reference 
body central force field and which may be expressed in closed form in 
terms of its true or eccentric anomaly. 

The above "two-body" dynamic model may be extended to include 
perturbational accelerations acting on the vehicle, in which case 
Equation (1) becomes: 



(?) 


where 


?\ is the summation of all the accelerations arising from 
planetary,lunar and solar attraction, 

?2 an acceleration arising from the oblateness of the 
reference body, 

— ► 

P 3 is an acceleration arising from drag as derived from 
an appropriate model to be described later, and 

— ► 

P 4 is an acceleration due to direct solar radiation upon 
the vehicle neglecting reflected sunlight from the 
reference bodies or planets. 

“► -► - ► 

The perturbational accelerations P^, P2, P3, and P 4 are obtained 
from specific dynamical models whose details are described below, 

PERTURBATIONS DUE TO PLANETARY ATTRACTIONS 

The general expression for the perturbational acceleration, Pj, 
of a space vehicle due to the gravitational influence of the Sun, Moon, 
and other planets (excluding the reference body) is given by: 


j 



(3) 


where 


UJ - GMj 

G is the gravitational constant. 


is the tnas 9 of the jth body. 
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R v j is the position of the vehicle with respect to the 
jth body, (Figure 1) 

—► 

Rrj is the position of the reference body with respect 
to the jth body, (Figure 1) 

AR is the position of the vehicle with respect to the 
reference body (Figure 1)* 


j th body 



Figure 1. Vectors for Planetary Attractions 


If the two terms of the bracketed expression in Equation (3) are 
very nearly equal t a loss of accuracy due to round-off errors by 
machine computations will occur. Equation (3) can be written in a 
more convenient form due to Battin uZfj thereby eliminating this pro¬ 
blem. Using Equation (3), and Equation (4) below and rearranging terms 
one obtains Equation (5). 

AR « R v j - R r j (4) 


- r hj + ,Rvj 3 . 

p ‘ ■ J R^r t«rj (^T - 1) - 4R 1 

This latter form is used to achieve more computational accuracy. The 
actual programmed equations, which can readily be shown equivalent by 
substitution are: 

Pi ■ l t^rj U(U)) - AR] 

where 

f(U ) - + u 

1 + (1+U) ^ 


(5) 


( 6 ) 


(7) 


U 


2(R rj + AR) • & 

R rJ 2 


( 8 ) 
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The forces introduced by solar, lunar and planetary attractions 
in orbit prediction about the Earth are usually smaller by a factor 
of 10“ 7 than that of the primary attraction of the Earth’s gravita¬ 
tional pull. Planetary perturbations (including solar and lunar per¬ 
turbations) do have a significant effect on a vehicle trajectory over 
long time periods or for deep space probes. 

PERTURBATIONS DUE TO EARTH OBLATENESS 

The fact that the Earth is not a perfect sphere of uniform den¬ 
sity gives rise to perturbational accelerations due to Earth oblate¬ 
ness. Formulas are derived below which the program uses to include 
these perturbations in the prediction of near-Earth trajectories. In 
a coordinate system attached rigidly to the Earth # the Earth oblate¬ 
ness perturbation may be treated as the acceleration of a conservative 
force derivable from a potential. Thus: 

P 2 - - VU (9) 

The potential function is given by Equation (12) and the final form of 
the acceleration ? 2 is given by Equations (26)-(29) and Equation 
(32) below. 

The Geopotential Function 

The geopotential function can be derived from the basic property 
that the potential, U, satisfies Laplace’s Equation^ 5 PP ’ 


V 2 U - 0. 


( 10 ) 


In spherical coordinates. Equation (10) becomes 


V 2 U 



au> ^ a r 
+ u ( cos * 


3U ^ 9 (■ 1 3U-| 
a<t>' ax 'cos $ ax' 


o (11) 


Although the solution of Equation (11) is not elementary, the equa¬ 
tion may be solved by applying the method of separation of variables 
and using Legendre polynomials ub ’ p,Zf0 -L The solution of Equation 
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(12) 


(11) is given by: 



n-o 


m*o 


where 


V is the Earth’s gravitational parameter, 

R e is the Earth's mean equatorial radius , 

r,A,<{> are spherical coordinates (see Figure 2), 

P™ is the associated Legendre function (spherical harmonics) 
of the first kind of degree n and order m, and 

Cn # m and 

s n,m are numerical coefficients. 

The Legendre polynomials , P n (x), are defined by 



and the associated Legendre function P n ( x ) by 


Pn(x) - (1 - x 2 )? P n (x). 


dx m 


The spherical coordinates r t A, are defined with respect to 
the fundamental planes determined by the true equator, and the Green¬ 
wich meridian (see Figure 2). 


z 


Greenwich meridian 



vehicle 


+ y 


true equator 


x 


Figure 2. Spherical Coordinates 
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The fundamental term In the expression for U is given by 
m ■ n ■ 0. The terms in which m ■ 0 are called zonal harmonics. 
Inspection of Equation ( 12 ) shows that these terms vary only with 
latitude, and hence reflect deviations of the Earth 1 s potential from 
a sphere of uniform density that are symmetric around the spin axis 
(e.g., a pear shaped Earth potential can be modeled with zonal harmon¬ 
ics) . Evaluating the 2,0 term for U we have 

- ^ (“t) 2 P2 ("1" *) C 2f0 

and since 

P 2 (x) - \ (3x 2 - 1 ) 

Equation (13) becomes 

- r (“r) 2 2 ° Sln2 ♦ ‘ X) C 2»0 • 

The zeros of Equation (13) are at * 35?25. 



Figure 3. Zonal Harmonic for n = 2, m = 0 

The sectorial harmonics arise when m - n. The zeros of the 
sectorial terms lie along lines of longitude. The remaining combi¬ 
nations are tesseral or square in that these terms vanish both along 
a number (m - n) of parallels of latitude and a number ( 2 m) of meri¬ 
dians of longitude. Figure 4 provides an illustration of zonal, 


03) 
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sectorial, and tesseral harmonic variation* for a sample set 



Zonal 

Harmonics 


Sectoral 

Harmonics 


Tesseral 

Harmonics 


m 


n 


- 4 

- 0 


n * 3 
m ■ 3 


n 


- 8 


Figure 4. Variations 


In summary, Equation (12) describes the Earth*s gravitational 
potential at any point in space; the negative gradient of this poten¬ 
tial gives the corresponding acceleration. Note should be made, how¬ 
ever, that the acceleration obtained is in the Greenwich coordinate 
system, and hence must be transformed into the inertial system (the 
Base Date System) to be employed in trajectory computation. 

Computation of Acceleration Dua jro -garth Oblateness 

The components of gravitational acceleration are now expressed 
as the sum of terms which have the same general form for any m, n 
combination. First, it is necessary to develop expansions for cos mX 
and sin mX in terms of cos X and sin X, From DeMoivre*s theorem. 


cos mX 4* i sin mX - (cos X 4* i sin X) m 


Now, expanding by the binomial theorem 



For m an even number 
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cos 


a 

2 


mA + i sin mA « £ (“ ) (-l) k (cos X) m “ /k (sin X) 


. m-2k, 


2 k 


k-0 


'S.1 
2 1 


k, , N m-2k-l, . , N 2k+l 
2k+l' x ^ cos ( sln X ' * 


+ 1 \ l ( 9V4.1 ) ( _1 ) ( COS 

k-0 


and for m an odd number 


m-1 

2 


cos mA + i sin mA - £ ((-l) lc (cos X) n ” ZK (sin A) 


v m- 2 k. 


?V 


2 


k-o 




m , N k, ^ v m-2k-l, ^ , N ?k+1' 


+ 1 S I (2k+l^ (_1) (cos X)I 


k-0 


(sin X)' 


Equating real and imaginary parts above, we have 


M 


cos 


k-o 

M ’ 


mX - l (-l) k ( ^) (cos X) 1 " 2k (sin X) 


m-2k 


2k 


sin mX - £ (-1) ( 2 k+l^ co9 X ' ~ _1 (sin A) - 


k-o 


no 


no 


where (®) - 


m! 


l J SL ! (m-£) ! 


if m is even 


M 


m -1 


if m is odd j 


f — - 1 if m is even 


and N f 


m -1 


if m is odd* 


From the definition of spherical coordinates: 
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sin <f> 


cos <f> cos X - — } 
r 


cos 


$ sin X ■ ^ 


(16) 


where x, y t and z are the coordinates of the vehicle in the Green¬ 
wich coordinate system. 

Using Equation (16) in Equations (14) and (15), we have: 


cos m\ - 


-L 


M 


(r cos <}>) m k-o 


c , 1N kf mi, .m-2k, .2k 

1 (- 1 ) l 2 kJ (x) (y) 


07. 


M' 


sin mX ■» 


l (-» k ( 2 ” +1 )w 


m-2k-l / N 2k+1 

(y) 


(r cos <j>) k-o 

Substitution of Equations (17) and (18) into Equation (12) yields: 

u c f R e1 n c ^ (® ln ^ 

UfJ l - - 


n-0 


m-n (r cos <J>) 


where 


M 


G(x, y) - C n#m l (-l) k (^ k )(x) 
k-0 


m-2k. N 2k 

(y) 


OS) 


(120 


M’ 


+ Sn.m l (-D k ( 2 " +1 )(») 

k-0 


m-2k-l, . 2k+l 

(y) 


From the definition of the associated Legendre function, we have 


m 


P >> 


x) 

2 n n! 


- d ®* 1 


dx 


m+n 


[(T 2 - l) n ]. 


m 


Setting t « sin ♦ and dividing Equation (19) by cos 0, 


09) 


P% (sin ♦) 

cos m 
n 


2 n! dx 


[(T* - 1)"1 

m+n 


Setting U - £ £ U m>n , and using Equation (20) we have 

n»0 m»0 


(20) 
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. _ 1L 
•“ r 


73 ' 

JL. 

n 

i 

d nTfn 

^ r ^ 


_r m 2 n n!. 

, m+n 
dx 


(t 2 - l) n H(x t y) 


-UR, 


m+n+1 , 

r 2 n! 


i m.+n 


dx 


m+n 


(t 2 - l) n C,(x, y) 


The gravitational acceleration is computed by taking the negative 
gradient of the potential function. The general derivative in the 
gradient will be taken with respect to £» where £ takes on the 
values of x, y t and z. 

Defining 


A 5 

Am t n 


au 


rum 




and carrying out the indicated operations using 


£r _ £ 

95 r 


and 


95 35 LrJ r |.H r 2 J 


yields: 


U R e n 


A n,n 


2 „ 


OU.y) (,* - 1)“ 

dx 


_ g(x.l y) [ii _ iil .dSULlI ( T 2 _ D n 

r 13£ r^J , m+n+1 


<T 2 - 1) 


dx 
n 3G (x 


dx 


m+n 


y)1 

95 J 


where 


3C . ( . |i .Zl - Cn,m l [(“D k ( 2 m k ){(m-2k) x m ~ 2k ~ 1 y 2k || 


95 


k-o 


35 


, m-2k 2k-l 8yi, 

+ 2k x y «^}] 


M * 

r> . . k f m ■>[/ m—2k—2 2k+l 3x 

+ S n ,m I [(-1) ( 2 k+lM ( ®“ 2 k-l) x y — 


k-0 


95 


m-2k-l _2k ^ 
95' 


+ (2k+l) x*™“~ y— ii}]. 


( 21 ) 


( 22 ) 
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that 


where 


It can be ahown by applying the binomial theorem and induction 

4 « 2 - »“ l w> k «) <’> 2n ' 2k ’ 1 

dT k-0 

[n - 4] 


and 


“] indicates the integral part of n - Hence, 


d m+n 


[ 2 = 2 .] 


( 2 _ j\«* m y / j\k r n- j ( T ) n- 2k-m 

. m+n V ’ 1 ' l) l k' (n-2k-m) ! U ' 

dx k-0 


( X 2 - i) n - l (_D k ( n ) l&E&H (x) n-2k-m-l 

m+n+1 V X ' <• V ■ L '' '•W fn-2k-m-l'»! KT J 


f m-n-l, 
l 2 J 


._m+n+l 
dT k-o 


k J (n-2k-m-l)! 


( 23 ) 
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where t - sin <J> - * . Therefore, Equation (21) may be rewritten as 


[2JL] 


*n,m 


1 R<> n "T* (_i)k r*»i (2n-2k)! (• j zsn-2k-mj 

_n , nrfn+1 L l L 1 ' l k J (n-2k-m)! V 1 

2 nt r k-o 


[S^azi] 


[ teap L 4 G(x, y) - - 


34 


i i <-» k is) (fr 2k ’”' 1 ) 


k-o 




(25) 


By separately considering the case for n~m odd and for n-m even, 
it may be shown that Equation (25) reduces to the following expres¬ 


sion: 

X 


[id] 

u Re n 1 2 J 


.4 u *e“ r i, , vk _.(2n-2kl!_ r ^n-2k-ra-l, M f?t> 

An » m ‘ „n nrfn+1 l t <- 1} k! (n-k)! (n-2k-m)! U 1 { 


2 r 


k-0 


{((2n-2k+l) ^ - toT , ^ |£) 0(x , y) - f iSiSjjl) 
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where 


(27) 


3x il if ( - x 3^ „ rl if 5 - y rl if 5 - 

35 *0 if 5 t x* 35 " ‘0 if 5 ¥ y* 35 " ‘0 if 5 ¥ 


and 


M 


G(x, y) - 


Gn|tn I (“1) 
k-o 


m! , N m-2k, N 2k 

“ ..;■ (x) (y) 


2k! (m-2k)! 


M 1 


,m-2k-l . . 2k+l 


, c v / n \k m; _ . *m-2k-l 

+ S n,m L (2k+l)! (m-2k-l)! (x) (y) 

k-0 


(2R) 


M 


m/2 if m is even 


m-1 


M 1 - 


| - 1 if m is even 


if m is odd 


m-1 


if m is odd 


M 


22^- c n<m I K-i) 


m! 


k-0 


2k! (ra-2k)l 


{(m-2k)(x) 


m-2k-l / x 2k 3x 
(y> — 


+ 2k(x)"- 2k (y) 21 '- 1 ii(] 

35 


M’ 


Sn t m I [(~1) 


ml 


k-0 


(2k+l)! (m-2k-l)! 


// , w N m-2k-2 / x 2k+l 8x , /01 N m-2k-l / .2k 3yi , 

X {(ra-2k-l)(x) (y) rr + (2k+l)(x) (y) -ff ] 

dt, o t. 


(29) 


It is important to note that when z - 0, a special procedure 
must be used since the term (“)” 1 arises in Equation (26) for n-m 
even; and (^)°* when n-m is odd. 

When n-m is even. Equation (26) reduces to one term, when z = 0. 

n-m 

V (n+m+l)S G 3G 

2 ;* \*— 2 "’ " ' Hl ' 


n r n—m 

.5 _ ~ y R e . [. 2 (n+m) 1 J f(n+m+l) 

n,n 2 n r m+n+ l V' (S=a)! (n-SCftjJl 


( 30 ) 


Similarly, when n-m is odd. Equation (26) may be reduced 
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^n,m 

An,m 


0 for t * x, y 


v» R e 

0 n m+n+1 
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■letsai!_) 
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( 31 ) 


In summary, the program computes the acceleration of a vehicle 
due to the Earth in the Greenwich coordinate system by the formula: 


00 n 

p 2 - I l <A£, m t + j + An, m k) (32) 

n -0 m *0 
~t 

where i, j, k are unit vectors in the Greenwich coordinate system. 

In Equation (32) the range of indices (n,ra) are restricted by 
SPACE-A to the following limits: 

n m 

2 0 , 1 , 2 

3 0, 1, 2, 3 

A 0, 1, 2, 3, A 

5 0, 1, 2, 3 

6 0 , 1 , 2 

7 0,1 

8 0 

9 0 

10 0 

The fundamental term is given by the (0, 0) combination and 
would be the only term present if the Earth were of uniform density 
and a perfect sphere. Since the center of the coordinate system is 
taken to be at the center of mass, it can be shown that the (1, 0) 
and (1, 1) combinations vanish^* b . 

F 

In Equation (32), A^ ^ is given by ( 26 ) when z 7 * 0, and when 
z = 0, by Equation (30) or Equation (31) depending on whether n-m 
is even or odd. 
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PERTURBATION DUE TO DRAG 


Accurate simulation of an artificial satellite or other space 
vehicle trajectories requires consideration of vehicle deceleration 
resulting from atmospheric drag. A number of planets (e.g., Earth, 
Venus, Mars and Jupiter) have sufficiently dense atmospheres to re¬ 
tard the motion of a vehicle within their atmospheres; however, this 
discussion confines itself to the atmospheric model of the Earth. 

An analysis of drag must take into account the particular mis¬ 
sion of the vehicle, e.g., low eccentricity, orbit, reentry, or fly¬ 
by, since vehicle mission determines what portion of the atmosphere 
it is necessary to include in the model. 

The following discussion describes the general equations used 
for drag computation, some of the problems involved in simulating the 
Earth's atmosphere, and the effects of certain simplifying assumptions 

Drag Equations 

The program uses two different formulas for computing the vehi- 
cle deceleration, P 3 , due to drag. For a relatively dense atmos¬ 
phere where the assumption of continuum flow is valid the following 
well known equation is used: 

-v 1 CpS + 

p 3 - - Cl [j P —] V a v a 

where 

p is the atmospheric density at the vehicle, 

Cd is the drag coefficient of the vehicle (dimensionless), 

S is the effective surface area of the vehicle, 

m is the mass of the vehicle, 

V a is the vector of the velocity of the vehicle with 
respect to the local atmosphere, 

V a is the magnitude of V a , and 

Cj is a constant used to convert the above expression 
into the units used by the program. 

Equation (33) is used to compute the acceleration of drag in the 


(33) 
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lower atmosphere from 0 kilometers up to 120 kilometers, although the 
model for the lower atmosphere from which the values of p and c l) 
are obtained may be extended up to 210 kilometers with some loss of 
accuracy. The basic model used in the lower atmosphere is the l?, S» 
Standard Atmosphere 1962^ 3 and its details are discussed belov:. 

As the atmosphere becomes more diffuse* the mean free path length 
(average distance between impacts of air molecules) increases. At 
110 kilometers, mean free path length is roughly one meter and at 130 
kilometers may become as large as ten meters 1 When the mean free 
path length exceeds the dimensions of the vehicle, the assumption of 
continuum air flow is no longer applicable. In such a diffuse atmos¬ 
phere where all collisions are essentially two-body collisions, the 
air ilow is referred to as free molecular flow. 

Ketchum L ^ J has derived the following formula for the magnitude 
of drag deceleration in free molecular flow: 

\h\ “ f (1 + Cav f] V a 

where 

R the radius of the vehicle, 

X the mean free path, 

C a v the average velocity of particles in the medium. 

Ketchura is uncertain as to the validity of the (1 + 2R/X) term 
in Equation (34). In the high atmospheric region the assumption that 
X >> R is usually justified, except perhaps below 140 kilometers. 
Therefore, the program actually uses Equation (35) in computation of 
the drag in the high atmosphere. 

P3 “ — ^2 [g*P Cav ~] ^a 

where, 

S is the effective surface area of the vehicle, 
m is the mass of the vehicle, 
p is the atmospheric density at the vehicle, 

C av is the average velocity of particles in the medium, 


(34) 


(35) 
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+■ 

V a is the vector of the velocity of the vehicle with 
respect to the surrounding atmosphere* and 

Cp is a constant used to convert the expression into 
the units used by the program. 

Equation (35) is used to compute the deceleration of drag in the 
upper atmosphere from 100 kilometers up to 1*050 kilometers* although 
the model for the upper atmosphere may be extended up to 2*000 kilo¬ 
meters, The basic model used for the upper atmosphere is that due to 

r 9 n 

Harris & Priester L The values of p and Cav used in Equation 

(35) are derived from this model, the details of which are also dis¬ 
cussed below. 

Notice in both Equations (33) and (35) that the direction of the 
drag force is in a direction opposite to the velocity with respect to 
air. In addition* both formulations assume zero lift and assume that 
the angle of attack of the vehicle is zero* i.e., that the vehicle 
velocity relative to the air mass is in line with the vehicle longi¬ 
tudinal axis. 

It is readily seen that the formula for drag in the upper atmos¬ 
phere* Equation (35),differs from the equation of drag in the lower 
atmosphere* Equation (33), Furthermore* in the region between 100 
kilometers and 120 kilometers there are disagreements between the two 
models. For example* the value of density predicted by the low atmos¬ 
pheric model (U, S, Standard Atmosphere 1962 ) at ]20 Vilorrters is 
34% higher than that predicted by the high atmospheric model (Varris- 
Priester) at the same height. 

The present program achieves a compromise between those two 
models by treating the region from 100 to 120 kiloreters as a tradi¬ 
tion region in which a weighted average is taken between the drar 
values computed by the two rethods and gradually slides the weight 
from unity for free molecular flow and zero for continuum flow at 120 
kilometers to unity for continuum flow and zero for free molecular 
flow at 100 kilometers. 
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because of the uncertainties in the atmospheric models and be¬ 
cause of the approximations made in the analysis the computation of 
drug deceleration is probably accurate to + 5% in the lower atmosphere 
and is less accurate in the upper atmosphere. 

Variables Used in Drag Computation 

Vehicle Mass 

In the most general case, the vehicle mass in the drag equations 
should be considered as variable with time. In the orbiting case or 
the fly-by case, a step change in mass representing the separation o; 
a landing craft is conceivable. A long-term steady-state mass flow 
rate, however, would probably be small. 

For the reentry case, if the reentry vehicle is of the heat-sim. 
type, the mass would be constant. For an ablative nose cone ^i.e. t 
one which loses mass when moving at high speeds due to friction) the 
mass flow rate is a function of the drag. For ballistic missile ap¬ 
plications, this mass change is usually ignored. In any event, such 
changes in mass represent a small error in the location of the impact 
point. Therefore, the program treats the mass as an input constant. 

Surface Area 

The effective surface area term S in the drag equation is not 
simply the cross-sectional area of the vehicle. The vehicle, in 
passing through the air, produces a shock wave which skirts the mis¬ 
sile thereby placing the effective cross-sectional area at a point 
somewhat close to the nose. Since the shock wave changes with air 
speed, so does the effective cross-sectional area. In practice, S 
is made constant and any variation with speed is included in the co¬ 
efficient of drag. 

Velocity with Respect to Air 

The velocity of the vehicle is available in an inertial coordi¬ 
nate system (Base Date System). 

—► 

The vehicle velocity with respect to the moving air mass, V a , 
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in the same coordinate system! is obtained by subtracting the velocity 
of the air mass from the vehicle velocity* A good first approximation 
to the velocity of the air mass is obtained by assuming the air mass 
to be rigidly attached to the rotating planet* 

From these considerations we have: 

v a - R - n' x r 

where 

V a is the velocity with respect to surrounding atmosphere* 

R is the position vector of the vehicle in the inertial 

frame t 

• w 

—► 

R is the velocity vector of the vehicle in the inertial 
frame* and 

-+ 

ft is the vector of angular rotation of Earth expressed 
in the inertial frame* 

A better approximation could be obtained by including the effects 
of wind velocity* The purely local effects have to be neglected* but 
the long-term horizontal effects are known as a function both of posi¬ 
tion on the Earth's surface and of altitude. The effects of the wind 
velocity's direction (independent of altitude but dependent on lati¬ 
tude and longitude) and magnitude (strongly dependent on altitude* 
less strongly on latitude* and least on longitude) would have to be 
included. The error made by neglecting Earth winds is about 1,500 
feet at impact for a typical ICBM mission. It should be noted that 
winds are of importance only in the Earth's lower atmosphere* mainly 
for the reentry case. In the present program the effect of winds is 
neglected. 

Drag Coefficient 

The drag coefficient* Cp* is sometimes considered to be a con- 

CdS 

stant. Very often the ballistic coefficient* B - * is used in 

* m 

analysis in place of the drag coefficient* mass* and effective sur¬ 
face area. A much more accurate representation for Cd is obtained 
by considering it to be a function of Mach number. Mach number is 


( 36 ) 
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defined to be: 



where, 

M Is the Mach number, 

|V a | Is the speed with respect to the surrounding air, and 
c is the speed of sound in the surrounding air. 

The speed of sound is a function of altitude obtained from the 
low atmospheric model (see below). It should be noted that as altitude 
increases, the atmosphere becomes rarified to the point that the speed 
of sound loses its physical significance. 

In the program Cd is tabulated for about 40 different Mach 
numbers. These numbers are denser for speeds below Mach 2 than these 
above and are very dense in the region around Mach 1, For intermedi¬ 
ate values of Mach number, linear interpolation is uaed. 

Inadequate knowledge of the drag coefficient is one of the major 
sources of inaccuracy in the simulation of drag. Since drag coeffi¬ 
cient is a function of Mach number, drag coefficient data has been 
obtained by wind tunnel measurements made at a range of Mach numbers. 

Air Density 

In the region below 120 kilometers the air density, p, is a 
function of altitude obtained from the low atmospheric model and is 
computed from a stored table (see below). 

In the high atmosphere p is obtained from the upper atmospheric 
model of Harris-Priester and is considered to be a function of alti¬ 
tude, local solar time, solar flux, and latitude. Tables are provided 
in the program for its computation (see below). 

Mean Particle Velocity 

The mean particle velocity, C aVt ie of importance only when 
the vehicle is in the upper atmosphere where the assumption of free 
molecular flow is valid. From Equations I,3,4-l and I,2,6-1 of 
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( 38 ) 


Reference [7], C is given by 



where 

T is the absolute temperature of particles ( # K), 
m is the mean molecular weight of medium (gra), 
k is a constant of proportionality (.IAS). 

The values of T and m are given in the Harris-Priester model of 
the upper atmosphere. C av can, therefore, be obtained directly from 
stored tables. 

Lower Atmospheric Model 

Drag in the lower atmosphere (below 120 kilometers) can be large 
and a vehicle entering this region will usually be slowed down suffi¬ 
ciently to be quickly captured by the Earth. Thus, the lower atmos¬ 
phere is primarily of concern in the reentry case. 

Data for an average model have been well established for the 
lower atmosphere. There are five sources for these data: U. S, 
Standard Atmosphere 1962 : COSPAR International Reference Atmosphere 
(CIRA), 1961; COESA Table for Tropical Latitudes, 1962; ARDC Model 
Atmosphere 1956, 1959. Table I shows the density deviation (in 
percent), as a function of altitude, of each of the others from the 
U. S. Standard Atmosphere values. From the table, it is evident that, 
except for the COESA tables, there is good agreement between the vari¬ 
ous tables at low altitudes. Note that the U. S. Standard Atmosphere 
and CIRA tables are in excellent agreement all the way to 120 km. 
(400,000 feet). 

The lower atmosphere is characterized by seasonal, diurnal, and 
latitude variations; however, none of these is sufficiently well 
documented. The only effect of omitting them is that the impact point 
of a re-entering body would be slightly different. It was estimated 
in 1958 that the standard deviation for a heat-sink type nose cone 
used in the ICBM application is only about 0.5 nm. 
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Table I 


Comparison of Sources of Density Data 


A1 

Lituce 

U. S. Standard 
Atmosphere 
Density 

Percent 

Devintir 

f ror ID : 

- 1 

:. ri **( 



V alues 
(Reference) 
slugs/ft 3 

ARDC 

1956 

ARDC 

1959 

n T d a 

cors.A 

1962 

km. 

ft. 

1 r\A 

1961 

0 

0 

2.38” 3 

0 

0 

0.55 

-4.77 

3.0 

10,000 

1.7b" 3 

0 

0 

-0.91 

-5.32 

5.3 

18,000 

1.36“ 3 

0.04 

0 

1.85 

-1.67 

IC.l 

33,000 

7.97“ 4 

0.05 

0 

1.68 

1.92 

14. b 

48,000 

4.00" 4 

0.09 

0 

2.36 

15.5 

20.4 

67,000 

l.bl“ 4 

3.28 

0.16 

0.48 

6.80 

29.0 

95,000 

4.20” 5 

0.59 

-2.36 

0.10 

0.46 

32.5 

110,000 

2.07" 5 

-3.13 

-3.13 

0.68 

2.53 

-.8.8 

lb0,000 

2.32 -6 

4.77 

4.77 

0.77 

8.92 

b7.1 

220,000 

2.50" 7 

15.0 

15.5 

1.30 

8.10 

SI.4 

300,000 

4.62 -9 

31.2 

-10.8 

0.11 

- - 

| 121,9 

400,000 

3.b2“ 11 

81.5 

-35.0 

1.17 

"" J 


Three tables with 32 values each are stored in the program for 
the lover atmospheric model. The first contains 32 values of alti¬ 
tudes in kilometers| the second contains 32 values of the logarithm 
(base 10) of the density, p, in gm/km 3 , and the third contains 32 
values of the speed of sound, c, in km/sec. The nth entry in the 
altitude table corresponds to the nth entry in the table for logjoP 
and c. Intermediate values are obtained by linear interpolation. 
These three tables are obtained from the U, S, Standard Atmos¬ 
phere 1962 , in which densities and speed of sound at altitudes below 
90 km are listed. Above 90 km the speed of sound was calculated from 
temperature and mean molecular weight data which are directly avail¬ 
able. 
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The value of altitude (in Earth radii) above an ellipsoidal 
Earth is obtained using a formula found in Baker^^. 

h - r - 1 + f sin 2 <J> + -y (7 - sin 2 2 <j> (39) 

where 

h is the height of the vehicle (ER), 

r is the distance of the vehicle from the geocenter (ER), 

4 > is the geocentric latitude of the vehicle, 
f is the flattening constant of the Earth « Vqr ' V • 


For programming purposes Equation (39) is put into the more convenient 


form: 


h - R - 1 + z l 


m ( 2 iw x 2 + y 2 ^ 
[j * [x ~ ~2^ 298.3 R 2 j 


298.3 R 2 


(40) 


where , 

h is the altitude of the vehicle (ER), 

x,y,z are the position coordinates of the vehicle in the 
Greenwich coordinate system (ER), 


2 is the distance of the vehicle from 
the geocenter (ER)• 

The position coordinates of the vehicle are available from the pro¬ 
gram and h is multiplied by 6378.165 to convert its units into 
kilometers. 

After h is obtained, p (p ■ 10^°® loP ) and c are obtained 
from the low atmosphere tables. Then the air velocity, V a and the 
drag coefficient, Cp, are computed as described above. Finally, 

P 3 is computed according to Equation (33). 

Upper Atmospheric Model 

Models of the Earth’s upper atmosphere (above 100 kilometers) 
must take into account solar activity. There is evidence that solar 
activity occurs cyclically at periods of 27 days, 6 months, 1 year 
and 11 years. 


R - /y 


+ + z 
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Theoretical models do not exist for the 27-day, 6-month, and 1- 
year cycles. Diurnal variations, if any, of the models for these 
cycles are not known. Investigation of the 11-year cycle (corres¬ 
ponding to the sunspot period) in solar flux led to the Harris- 

I 9 111 

Priester model of the upper atmosphere. This model L * J has 
diurnal and solar flux variations. 

r 91 

The Harris-Priester model J lists the density, absolute tem¬ 
perature, and mean molecular weight of the atmosphere as a function 
of altitude, solar flux and local solar time. The mean particle 
velocity can be found by use of Equation (38). Note that at the 
North and South Poles local solar time is undefined. 

The upper atmosphere has a delaying effect on solar radiation. 

It takes several hours for the Sun’s heat to pass through the atmos¬ 
phere and reach the Earth's surface. The Harris-Priester model is 
based on densities computed at the Earth's equator. Intuitively, it 
is expected that it will take longer for the solar flux to reach the 
poles as opposed to the equator. Therefore, it is considered that 
there is an effective variation of solar flux with latitude. This 
variation is implemented in the program by applying the Harris-Priester 
model at the equator and a stored table of "twilight" densities at 
the poles . The cosine of the latitude of the vehicle is used as a 
weighting factor to interpolate between the two sets of data. 

The Harris-Priester upper atmospheric model has been incorporated 
in the program by means of a table look-up procedure.. Two tables are 
used in the program. 

The first table lists the logarithm to the base 10 of the density 
in gm/km 3 times the mean particle velocity in km/sec (logjQP C a v) a t 
the equator as functions of altitude (km), solar flux (watts/m 2 - Hz 
at a wavelength of 10.7 cm) and local solar time (hrs). These values 
have been tabulated for 16 values of altitude, 4 values of solar flux 
and 13 values of local solar time. Hence, this table has 832 entries. 
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The second table lists logic pC 8 v * n t ^ ie 8aJne units at "twi¬ 
light" for the modeling of the polar region as a function of altitude 
and solar flux. This table has 16 entries for altitude* A entries 
for solar flux, or 64 entries. 

For intermediate values of the input variables (altitude, solar 
flux, local solar time) a linear interpolation is used to obtain the 
output, logic PCav« This method gives fairly accurate results since 
oCav is nearly exponential. The value of solar flux is determined by 
input data. The value of altitude is computed according to Equation 
(AO), The value of local solar time is computed from the x and y 
coordinates of the vehicle and the Sun in the inertial (Base Date) 
coordinate system (see Figure 5), 



Figure 5. View of x-y Plane from Above 

Thus we haves 

local solar time "{(02 - 0i + *rr) + 720°} (mod 360°) 


where 

01 
0 2 


tan" 1 (“) 
tan" 1 {^] 


where tt ^ 0 i ^ ~ it 
where * i ®2 i " 11 
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Xy t y v are the vehicle coordinates in the inertial system, 
x s* ys are the Sun’s coordinates in the inertial system* 

Local solar time given above is then divided by 15 to convert its units 
to hours* 

After a value of logiQ pC av has been obtained from both the 
equatorial table and the ’’twilight 11 polar table an interpolation is 
done to obtain the final value of pC a v using the cosine of latitude 
as a weighting function* Hence: 

L * Lp + cos $ (Le “ Lp) (J ''' 


X 2 + v 2 

2 + y2 + z 2 

pCav " 10 

where 

L is the final value of logiQ pC a y, 

Lp is the value of logiQ pCav obtained 
from the twilight table, 

Le is the value of logiQ pC a v obtained 
from the equatorial table, 

x,y,z are the coordinates of the vehicle in the inertial 
(Base Date) system and $ is the latitude of the 
vehicle* 

Once pC a v has been found, the vehicle area and mass (S and m, 
standard program inputs) as well as air velocity V a (see Equation 
(36)), are obtained and the value of drag deceleration is finally 
computed according to Equation (35). 

PERTURBATION DUE TO DIRECT SOLAR RADIATION 

Solar radiation exerts a pressure on the intercepting surface of 
a vehicle. Orbiting planetary vehicles, having a large area to mass 
ratio are subject to noticeable perturbations due to solar radiation 
pressure. In fact, for orbits above 500 miles the solar radiation 



(43) 

(44) 
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perturbation is usually more significant than that due to drag L 
The vehicle acceleration due to solar radiation pressure depends on 
the area to mass ratio of the vehicle as well as the intensity of the 
Sun’s incident power at the vehicle and the fraction of solar illu¬ 
mination on the vehicle* The illumination factor is considered in 
three distinct regions in the following analysis: full sunlight, 
penumbral illumination, and no illumination (i.e., umbral region). 

In computing the solar radiation perturbation, P 4 , this analy¬ 
sis neglects the dispursive effects of planetary atmospheres which 
complicate the geometry of the umbra and penumbra* The analysis also 
neglects the effects of reflected sunlight from the reference body 
or any other planet. 


Acceleration Pue to Solar Radiation Pressure 


An expression for the acceleration due to solar radiation pres- 

[121 

sure given in Wolverton J is: 


where. 



p is the illumination factor, 

q is the reflectivity coefficient, 

A is the area of vehicle pertaining to solar 
radiation pressure, 

m is the mass of vehicle, 

Lq* 3*86 x 10 26 watts, the total power output 
of the Sun (+ 3%), 

c is the speed of light, 

R sv is the vector from the Sun to the vehicle. 

The SPACE program assumes a reflectivity coefficient equal to 
unity. This assumption may be altered by changing the area. A, by 
an appropriate factor. Simplifying Equation (45) and putting the 
above variables into units used by the program, one obtains: 



(45) 


( 46 ) 
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where, 

p is the illumination factor, 0 < p ^ 1, 

A is the area pertaining to radiation pressure (ft 2 ), 
m is the mass pertaining to radiation pressure (lb-masses), 
Rsv is the vector from the Sun to the vehicle (ER), 

Cp - 1.04819 x 10 3 and a constant used to convert 


the expression into units used by the program and to ac- 
LO 

count for the —“ factor of Equation (45). 

Equation (46) is the basic equation used by the program for com¬ 
puting the solar radiation perturbation acceleration. 

The illumination factor, p, is obtained from the relative geo¬ 
metry of the vehicle with respect to the Sun, the Mocn and the planet'", 
p - 1 in full sunlight, 
p - 0 in the umbral region, 

0 < p < 1 in the penumbral region. 

It is possible that a vehicle may lie within the penumbral region 
of two bodies at the same time, e.g., the Earth and the Moor v In 
this case, only the penumbral illumination factor due to the reference 
body is computed# Errors introduced by this assumption are extremely 
small; first, because for most vehicles of interest the solar radia¬ 
tion perturbation is itself small (usually less than 10 5 times the 
acceleration of the reference body except for low density balloons); 
second, because only a short time is spent in the penumbral region of 
the reference body by an orbiting vehicle; and third, because the 
incidence of simultaneous penumbral obscuration is rare. 

Since the geometry used in calculating the illumination factor 
is the same as that used for eclipse information, the portion of the 
program which computes the solar radiation pressure perturbation also 
computes the times at which the vehicle enters or leaves the umbra or 
penumbra of the celestial bodies. The geometry used in calculating 
the vehicle illumination factor, p, is described below. 
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Vehicle Illumination Factor 


Figure 6 illustrates the geometry of vehicle illumination, neg¬ 
lecting the effects of atmospheric refraction* By similar triangles 
it is seen that the height of the umbral cone (h u ) and penumbral 
cone (hp) are respectively given by: 


where. 



R S p is the distance from the center of the Sun to the 
center of the reference body, 

R s is the radius of the Sun, and 

Rp is the radius of the reference body. 

Next, criteria are developed to see in which region the vehicle 
lies, i.e., full sunlight, umbra, or penumbra. First, consider 
Figures 7 and 8 and the definitions and relations below. 

R S p is the vector from the center of the Sun 
to the center of the reference body, 

R is the vector from the center of the 
reference body to the vehicle. 


? 


h u 


^8P 

Rsp 


( 47 ) 


(48) 


(49) 


Q - - h p 

* K sp 


From Figure 7 we get the following: 


(50) 


c ° 8 a ■ A - © 2 


cos a 


(R - P) » ^8P 

|£ - ?l Rsp 


(51) 

(52) 
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Figure 6. Vehicle Illumination 
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Figure 7. Umbral Region Geometry 
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Figure 8. Penumbral Region Geometry 
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From Figure 8 we get; 


cos 



( 53 ; 


C». b - ( * - 1* 8 r <«; 

|£ - Q! Rap 

Nov, if the scalar product then the vehicle 

lies on the side of the reference body away from the Sun. In this 
case, if cos a > 0 and if a A, the vehicle or satellite lies 
m the unbia and p is set equal to zero. Another test is: 

if | cos a | ^ | cos A | (*>b 

then p “ 0, 


i*e., the vehicle lies in the umbra. 

If the vehicle does not lie in the umbral region one can test to 
see if it lies in the penumbra. Thus, if cos 6 > 0 and 0 B the 
vehicle lies in the penumbra, or equivalently: 

if |cos 0| >, |cos B| (56) 

then 0 < p < 1 

£i„e., the vehicle lies in the penumbra). 

If the vehicle does not lie in either the umbral or penumbral 
region, then the reference body does not obscure the Sun's rays. A 
check can then be made to see whether any other celestial body blocks 
or partially blocks the solar radiation; and if not, the illumination 
factor is set equal to one, p ■ 1. 


Penumbral Illumination Factor 

If the tests of Equations (55) and (56) above indicate that the 
vehicle lies in the penumbral region, then the illumination factor, 
p, is computed according to the formula: 


P 


„ Ae * 
6 c* 


(57. 
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where Aex Is the angular area subtended by the exposed portion of 
the solar cap at the vehicle position, and 

65 is the total angular area of the solar cap at the 
vehicle position. 

The general approach used by the SPACE program for computing the 
penumbral illumination factor is somewhat more detailed than Chat 
given in most references, It consists in projecting the solar disc 
(or cap) and the cap of the obscuring planet (or moon) on to a great 
imaginary sphere whose center is at the vehicle position. The rela¬ 
tive angular areas of the caps are computed and the illumination fac¬ 
tor is given according to Equation (57), While the theoretical de¬ 
velopment is somewhat involved, it leads to simple closed form alge¬ 
braic expressions convenient for use on a computer. 

First, consider Figure 9 which shows a sphere of radius R, re¬ 
presenting the Sun or a planet, and a point P at a distance l 4* R 
from the center of the sphere. The apparent angular area of the 
sphere as seen from point P is equivalent to the angular area of 
the sphere’s cap projected onto a projection sphere of radius a. 
Therefore the following is trues 


angular area 0 - 2 tt [1 - cos 


- 2tt 


. h a + 2R) 

L * + R . 


where a is the radius of the projection sphere, and 
H is the height of the projected cap. 

The total angular area of a planet (0p) or of the Sun ( 6 g) at 
the position of a vehicle can be obtained by substituting the appro¬ 
priate values of l and R into Equations (58) and (59) and simpli¬ 
fying. Thus, we have: 


(58) 

(59) 


9p *» 2ir 


/h(h + 2Rp) 


1 - 


h + R P 


- 2n 


H’ 


(60) 
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Figure 9. Geometry Used to Measure the Angular Area 
of the Solar or Planetary Cap 
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I 


(61) 



where 

0 p is the angular area of the planet (or moon) at the vehicle 
position, 

h is the height of the vehicle above the planet’s surface, 

Rp is the radius of the planet, 

0S is the angular area of the Sun at the vehicle position, 

R S v Is the distance from the center of the Sun to the vehicle, 

R s is the radius of the Sun, 

H is the height of the solar cap, 

H f is the height of the planetary cap, 

a is the radius of the imaginary sphere of projection* 

Equations (60) and (61) establish the relative size of the solar and 
planetary caps as seen from the vehicle* 

Next, it is desired to determine the angular area of the exposed 
solar cap Agx» To do this, consider Figure 10. In this figure the 
vehicle is positioned at the center of a great imaginary projection 
sphere of arbitrary radius, a* The relative positions of the Sun 
and the obscuring body or planet are shown as well as the projections 
of the solar and planetary caps onto the great sphere. The equator 
of the great sphere is constructed to be coplanar with the vehicle, 
the center of the Sun, and the center of the planet. A great circle 
is also constructed perpendicular to the equator and passing through 
the intersection of the caps; this circle will be used to define one 
of the limits of integration in the computation of the angular area 
of the two lunes which are formed by the great circle. 

For the calculations to follow, three coordinate systems are de¬ 
fined and illustrated in Figures 10 and 11. The coordinate systems 
employed are: 
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and Planetary Caps 


y* y 1 * y H 



Figure 11 . Definition of Coordinate System 









Cl) A system (x,y,z) for ths solar cap, 

C2) A system (x , ,jr , ,z') for the planetary cap. 

(3) A system (x",y",z") for the great circle. 

All three coordinate systems are orthogonal and have. In common, the 
same y axis, The relations between the coordinates are given In 
Equations (62) and (63 )t 


x' 


cos 0c 

sin 0C 

0 


X 

z' 

- 

sin ©c 

cos 0c 

0 


Z 

_y’_ 


0 

0 

1 


y 

"x”~ 


sin 0 q 

cos 0 g 

o' 


X 

z” 

- 

- cos 0c 

sin 0 q 

0 


z 

y" 


0 

0 

1 


y 


where 

6c Is the angle from the center of the solar cap (z-axls) to 
the center of the planetary cap (z'-axis), 

0G is the angle from the center of the solar cap (z-axis) to 
the great circle (x"-axis). 

By inspection of Figure 10 one finds that: 


(62) 


(63) 


cos 6c 




R sv 


(64) 


where 

0 < 6c \ 

Next find the points of intersection of the planetary and solar 
caps x ■ Xp f y ■ yp f and z - zp. Referring to Figures 9 and 10 
one sees that the equation of the small circle of the solar caps 
(having a cap height of H) is! 
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x^ 4* m 

z - a - H (65) 

. \ x 2 + y 2 - H (2a - H) (66) 

Similarly! the equation of the small circle of the planetary cap 
having height H f is: 

(x') 2 + (y*) 2 + (z') 2 - a 2 

z f - a - H f (67) 

x »2 + y»2 m (2a - H f ) (61 

Using Equation (62) we can write Equations (67) and (68) in terms of 

the (x,y,z) coordinates. 

z 1 * x cos 6q + z sin @c “ a - H* (69) 

x' 2 + y'2 - (x cos 0 C - Z sin 0 C ) 2 + y - H'(2a - H') (70) 


From Equations (65) and (69) we find: 

(a - H 1 ) - (a - H) cos 

y ■ y m MPa 

X X P sin 6 C 

Using this result and Equation (65), Equation (70) gives: 


(71) 


2 2 

y - yp “ 


H*(2a - H') sin 2 0 c ~ [(a - H*) cos 6 C ~ (a - H)) 

Bin 2 0£ 


(72) 


And, of course, Equation (65) can be rewritten: 


z m Zp » a - H 

Equation (72) requires some interpretation. Figure 10 illustrates a 
case where the solar and planetary caps intersect. It is possible 
that the caps are tangent or that the planetary cap lies within the 


( 73 ) 
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solar cap. To determine whether an intersection exists, the value 
of yp 2 of Equation (72) may be used as a discriminant. 


(1) 

If 

yp 2 

> 0, 

two intersections exist. 


(2) 

If 

yp 2 

- o, 

the caps are tangent. 


(3) 

If 

yp 2 

< 0, 

the caps do not intersect and 
exposed solar angular area is 

Aex “ e S - 0p. 

the 

given by 


For the great circle passing through the intersection of the 
caps z" - 0, z - zp, x - xp. Hence, Equation (63) gives: 

z" - - x cos 0 G + z sin 0 q - 0 (74) 

x X P 

° r 7 " “ " tan e G (75) 

thus from Equations (71) and (73) 


(a - H 1 ) “ (a - H) cos 0^ 
tan 6 G “ (a - H) sin 6 C 


In a similar manner we may also find the equation for 
the angle between the center of the planetary cap and the great 
circle^ 

tan 6 G ' - p- 

everywhere on the great circle. 

Using Equation (62) evaluated at x ■ xp and z ■ zp we 
finally obtain: 

(a - H') cos 0 q - (a - H) 
tan (a - H*) sin 0 q 


In computing Agy* there are five cases of interest and these 
cases are illustrated in Figure 12. If yp 2 < 0 one obtains case 
V where Agx ■ 0g - 6p. If yp 2 > 0, however, the caps intersect; 
and before AgX can be determined, one must calculate A s and 


(76) 


(77) 


(78) 
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Figure 13, Area of Spherical Lune, 
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Ap• where A s is the angular area of the smaller lune formed by the 
solar cap and the great circle, and Ap is the angular area of the 
smaller lune formed by the planetary cap and the great circle. 

Figure 13 illustrates the geometry and the variables used to 
compute A s , The method of computing Ap is identical except that 
all computations are done in the (x f , y f , z 1 ) coordinate system 
instead of the (x, y, z) coordinate system. 

The area of the surface ABC of Figure 13 is given by: 


area of ABC - 


*0 e b («) 

ad4> [ a sin 6 d6, 

0 M*) 



To get the angular area, A s , one must double the area of ABC and 
divide the result by a 2 . Thus, 


(79. 


As 


2 


* 0 6 b(*) 

d$ sin 6 dG 


0 $ a (4>) 


(80) 


or 


where 


A s 



[0 a (4 ) )] ~ cos [ e b(4 > )l 


d<J> 


0 


(81) 


4) is the angular displacement in the x-y plane, 
positive counterclockwise, 

0 a (4>) is the angular displacement of side BC (an 
arc of the great circle) from the z-axis, 
positive counterclockwise, 

0g(<$>) is the angular displacement of side AC (an 
arc of the small circle of the solar cap) 
from the z-axis, positive counterclockwise. 

Before Equation (81) can be evaluated, expressions for cos [0 a (4>)] 
and cos [0g(4>)] must be found. The equation of the email circle of 
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the solar cap is: 


x 2 + y 2 - H (2a - H). 


(32) 


Using the following transformation: 

x ■ a sin 0 cos <J> 
y ■ a sin 0 sin 4> 
z ■ a cos 0 

Equation (82) becomes: 

a 2 (sin 2 0 cos 2 $ + sin 2 0 sin 2 $) - II (2a - H) 


(83) 


or 


( 1 - cos 2 e) - H (2a - H) 


(84) 


For the small circle of the solar cap, 0 ■ 0^, and Equation (84) 
becomes 



( 85 ) 


,2 - a 2 


( 86 ) 


The equation of the great circle is 

(x") 2 + y J 

and using Equation (53) for x" we obtain: 

(x sin 0 q + 2 cos 0 g) 2 + y 2 “ a 2 . (87) 

From Equation (75) we note that everywhere on the great circle 

2 tan 0 g ■ x, (88) 

Hence, after substitution of Equation (88) into Equation (87) we have 

(2 tan 0 g sin 0 g + 2 cos 0 g) 2 + y 2 ■ a 2 . (80) 

Inserting the polar coordinate relations from Equation (83), yields after 
some reduction: 

cos 2 0 (1 - cos 2 0q sin 2 $) ■ cos 2 0g cos 2 <t> (90) 

or, since 0 - 0 a ($) along the great circle 
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COS 0r COS 


cos 0_ (<f>) 


J 1 - cos 2 0q sin 2 # 


Now we can evaluate A g of Equation (61) using the value of cos 0 a 
and cos 0^ above, 

4> n 

0 f cos Q r cos 

- a - H 


A s ■ 2 


r- 


cos 2 0g sin 2 # 


d# 


(92) 


Rearranging we have: 

$0 


A s - 2 


f COS * dt 2 f (1 -£) d, 

' J sec 2 0Q - sin^4> o 


Letting sin $ ■ x # c 2 ■ sec 2 ^ 


sin 4> 0 


Ag ■ 2 


dx 




- 2 fl -*) 


$0 


■ 2 sin 


^sin 4> 0 

f 


sec 0 q 


- 2 (1 - 7) # ( 


- 2 sin -1 [sin # 0 cos 6 g1 * 2 # 0 ( 1 “ 7 ) 


(93) 


From Equation (82) and Figure 13 f it is evident that 


xp 


cos *0 “ / 2 2 " /H (2a - H) 


/ x p 2 + yp 


yp 


sin #o “ 


yp 


y X p 2 + yp 2 (2a - H) 


, yp > o 


(94) 


( 95 ) 
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therefore, 


<fo ■ sin 


yp 

Al (2a - H) 


0 1 *0 1 7 


(96) 


Using Equations (93), (95) f and (96) the value of A s is obtained. 
Notice that the above development was done under the restriction, 

0 _ *o 1 in which case the entire lune of intersection, A s , is in 

the first and fourth quadrants. If the intersection of the two caps 
occurs in the second quadrant (as in Case I of Figure 12) t then <f> 0 
is in the second quadrant and the entire lune of intersection lies in 
the second and third quadrants. Instead of recomputing this case, 
note that (see Figure 13) if the lune of intersection, A s , lies in 
the second and third quadrants one could rotate the x-y plane about 
the z-axis by 180°, resulting in the problem that has just been ana¬ 
lyzed. Therefore, Equations (93), (94), and (96) give the correct 
value of the angular area A s whether actually lies in the first 

or second quadrant. Notice, however, that if 0 ^ T then 

xp * 0 and cos > 0; but if y <J> q ^ it then xp < 0 and 

cos $o < 0. Thus, cos 4> 0 can be used tc discriminate between the 
two cases just mentioned, 

A completely parallel development is used in order to compute 
Ap. The major difference is that all calculations are done in the 
(x'y'z 1 ) coordinate system. The resulting equations are: 


A P 


2 sin -1 [sin 4>o’ cos Bq 1 ] - 2 $o' (l 



(97) 


cos o * 


Xp cos 00 — (a ■“ H) sin Gq 
v'h' (2a - II') 


(98) 


sin 4*o * 


yp 

■fa' (2a - H') 


( 99 ) 
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1 

$ o 1 * sin 1 

yp 

7T 

/li' (2a - H' ) 

o i <Jo i - 

1 


where the primed quantities are related to the planet and Of; 1 is 
given by Equation (78). 

Finally, note that many quantities are expressed in terms of a, 
the radius of the great sphere, and H or H* the depth of the solar 
or planetary cap. B) use of Equations (60) and (61) we have 


!H 

a 


2tt 


H* b P 
a * 2 tt 

Thus all quantities may be expressed in terns of ' s and 0p which 

are already known. This is most easily accomplished by setting n « 1 

0 s 

(since the radius oi the sphere is arbitrary) and letting 11 *= — 
b P t7r 

and H’ - — . 

The computation of the exposed angular area of the solar 

cap is shown for each of the five cases of Figure 12 in Table 11. Note 
that the discriminnts used to determine which case is to be evaluated 
are yp', cos $q, and cos 4 q*# 

Finally the penumbral illumination factor is given by Equation 

(57); 

a EX 

p = T7 ’ 


(ion) 

(in; 

(10j; 
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NUMERICAL INTEGRATION 


The equations of motion for a space vehicle are second-order 
differential equations which describe the accelerations arising i ror. 
the forces acting on the vehicle. Accounting for the primary gravi¬ 
tational field of the reference body and four types of perturbative 
accelerations, the equations of motion (from Equation (2)), 


„ . . _ + Pl + p 2 + Pj + ^ (2) 

If the perturbations are considered to be zero, the right hand 
oiue of Equation (2) reduces to one term and the vehicle will follow 
a Kcplerian orbit which may be described in closed form in terms of 
its true or eccentric anomaly. Usually, however, part or all of the 
perturbations are included and Equation (2) must he numerically inte¬ 
grated , 

There are two basic methods by which the integration of Equation 
(2) may be formulated, Encke’s .method and Cowell’s method. If Equa¬ 
tion (2) were to be numerically integrated in a straight-forward 
manner, the integration would be known as Cowell’s method. The sim¬ 
plicity of this method is offset by the large accelerations which must 
be integrated. As a consequence of the acceleration magnitudes, small 
time increments have to be used in the integration, and machine round¬ 
off error accumulates rapidly. Independent evaluations at many com¬ 
panies and universities have shown that Cowell’s method requires more 

t u , r n , [10, pp. 228-242] 

machine time than other perturbational schemes. See Baker J 

tor a further discussion of Cowell’s advantages and disadvantages. 

Despite its drawbacks, Cowell’s method is still widely used and 

is especially suited if the accelerations due to perturbations are 

as large or larger than the term due to the central force field (e.r., 

during reentry). The program contains both methods and the option is 

left to the user. 

Historically, Encke’s method is older than Cowell’s although 
the former is more sophisticated. Cowell’s method requires a modern 
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high-speed computer to be practical, whereas hmi.e's was developed 
for hand computation. In Lncke's method, it is assumed that the per¬ 
turbative accelerations, Pj, are small compared to the reference 
body acceleration. Consequently, when drag accelerations arc small, 
the solution of Equation (1) is a good approximation to the true orhit 
Under these conditions, it is only necessary to integrate the differ¬ 
ence between the accelerations on the two-body orbit and the total 
accelerations acting on the vehicle. The equations of motion then 
become second-order differential equations describing the acceleration 
differences. Let 

l m R - Rtb 

•¥ 

where R^b is the position of the vehicle in terms of the two-body 
orbit. Then, 



Equation (10A) is integrated to obtain £ and These quantities 

are then added to Rtb and RTBt respectively, to obtain the instan- 
taneous position (R) and velocity (R) of the vehicle. The quan¬ 
tity { is commonly referred to as the ,, Encke" term. 


Encke's Method 

Rewriting Equation (10A) we have. 


l 



rtb + 








The above equation could be integrated directly; however, the 
Bt r 

term ^ 1 - tt-} is not well suited for the numerical calculation 
3 , 

since R^g/R 3 is very close to 1, Hence we rewrite this term by 
setting 


(103) 


(104) 
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<1 - 7T 
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R • R * *f (R * RtB " R * Rtb^ "* * ^ TR 

^B 

R • (Rtb + “ r tb ( r TB + R ) f • (2 Rjb + c ) 
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3 L , 
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k tb [(i + q) ® + l] i + O')" 


expanding tbu nuneratrr 


R- 


R TB 


- 1 


3q 


1 - 4 - 


9 7 

3q + C 

— • 

(l+'O * 


In summary, LncKe's method computes the deviations from 
orbit by integrating 


f3q + 3q 2 + q 3 " - 
i- FT ) R TB “ *» 

L 1 + (3+n) 2 J 


V ► 

i Pi 


the nomine 1 


where 


(2R tb + S) 

R^B 


Determination of Kenlerian Orbit Vectors 

Lncke’s method requires the calculation of the position and velo- 

* • 

city vectors, Rtb and Rtb* respectively, of the nominal Keplerian 
orbit. The position and velocity at time t may be written in terms 
of tne position and velocity at an earlier time to as follows: 

• 

RtbCO “ f RtbUo) + g Rtb(cq) 

and 

• • 

Rtb( t) - l Rtb( to) + s RTB(t 0 ) 


where f and g are explicit functions of the differential eccentric 
anomaly of the Keplerian orbit. The program uses Herrick's method to 
determine f, t, g, and g, and then computes Rtb( 0 and R TB(t) 


(105) 

(10b) 


(107) 

(108) 
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from Equations (107) and (108), The equations for computing f, f, g, 
and g are given below; their derivations may be found in Battin c ^ 
or Baker*- 10 ]. 


where 


where 


U - X J [ 


f - 1 - C/r 0 

(109) 

r o - 1Rtb( t o )| 

(110) 

M - U 


8 “ /— 

/u 


. - y~\T s 

A r Q 

(111) 

g - 1 - C/A 

(112) 

/T (t-to) - ro X + d 0 C + C 0 U 

(113) 

r-L . JLi + _ JiL + i 

l 2! 4 !a 61a 7 FTa 7 * * * J 

(114) 

r J- _ JLi + _Si- + , 

3! 5 !a 7 ! a 2 9!a3 * * ,J 

(115) 

C v U 

S « X - “ 

a 

(116) 

r TB( c o) • ^BCto) 

d ° " J— 

V \i 

(117) 

A ■ Tq + d()S + coC 

(118) 

r 9 v ° , 

c o - -- l 

(119) 

i. _L _ li 

(120) 


a r 0 

a is the semi-major axis of the Keplerian orbit, 
p is the universal gravitational constant. 


and 
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v 0 • I^Ti(to)I• 


Solution of Equations (109) through (112) requires the determi- 
nation of X. Equation (113) la aolvad for X by writing 

F(X) • M - r 0 X - d B C - c 0 U - 0 

and 

F’(X) - - r 0 - d 0 C' - c 0 U’ • 0. 

Salacting an initial estimate of X to ba Xg t tha root of 
F(X) • 0 may be detarmlnad by a Nawton-Raphson itaratlve process, 
i . e • f 

x i+l • Xi - F(X i ) /y'(Xi) 

From Equations (114) and (115) 

C' - X - U/a 
U' - C. 

Slnca Equations (114) and (115) ara sariaa expansions in powars 
of —, a method must be found to limit tha alza of this tevm, For 
an ellipse, 

X - (E - E 0 ) /r 

where E, and Eo ara the eccentric anomalies at times t, and to 
respectively. Hence, by updating tha epoch to at frequent intervals, 

the difference E - E 0 is kept email. By similar reasoning, fra* 

X 2 

quant updating will keep small for hyperbolic orbits. 

CowjllJ^Method 

As described before, tha general equations of motion of a space 
vehicle are: . 

i - I h <i2i) 

i-i 

In Cowell's method, these equations are integrated, using numeri¬ 
cal techniques, to obtain the instantaneous position and velocity of 
tha vehicle. 

5t 



The program performs the integration using the same techniques as it 
does for the Encke method. The Runge-Kutta starting procedure pro¬ 
vides the initial data, and the Nordseick method is used as the long¬ 
term integration procedure. 

Integration Technique 

Equation (2) in Cowell's method or Equation (105) in Encke's 
method must be numerically integrated by the program. This process 
is divided into two stages, a starting procedure and a long-term pro¬ 
cedure. The long-term numerical integration procedure requires know¬ 
ledge of previous data points. Thus, the starting procedure is needed 
to provide the initial data points for the long-term procedure. 

The procedure chosen by Sperry-Rand for long-term integration was 
based on their experience with methods used in the ITEM and MINIVAR 
programs, and the results of comparative testing. The justification 
of their choice is presented below as it appears in the original docu¬ 
ment on the SPACE program^ -1 . 

"The long-term numerical integration procedure presently in use in 
the ITEM and MINIVAR programs is an Adams sixth-order predictor method 
(without corrector) for second-order differential equations. It was 
desired, however, to test a broader class of procedures before deciding 
on one for use as the long-term numerical integration procedure to be 
used in the program. Accordingly, a program was written to test pre¬ 
dictor, predictor-corrector, and iterated predictor-corrector, i.e., 
repeated application of correctors, techniques of various orders of 
approximation and with or without modifiers. 

,, The following results were obtained: 

(a) Modifiers were found to leave the error unchanged. There 
is, therefore, no reason to use them. 

(b) As time interval increases, there is more tendency for the 
solution to become unstable. Error increases with a large 
power of the time interval. 

(c) As the degree of the approximating polynomial increases, a 
decrease in stability is noted. Error decreases about 3:1 
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for unit increase in degree of the approximating polynomial. 

(d) Predictor-only methods (degree and time interval held fixed) 
are about 30:1 less accurate than predictor-corrector meth¬ 
ods! i.e.p one application of the corrector removes 97% of 
the error in the predictor. A second iteration of the cor¬ 
rector does not reduce the error, but does improve stabil¬ 
ity at the expense of a 2:1 increase in running time. The 
increase in running time is intolerable; therefore, a pre¬ 
dictor-corrector method will be used with no iteration of 
the corrector. 

(e) Either a fifth or sixth degree approximating polynomial 
yields a good compromise between accuracy and stability. 

"The final choice of a long-term integrating procedure involves 
considerations other than only accuracy and stability: 

(a) The ease of transforming the output of the starting proce¬ 
dure to the form of input starting data needed by the long¬ 
term procedure. 

(b) Whether the long-term procedure can easily accomodate a 
change in the time interval. 

(c) Whether the long-term procedure can easily interpolate 
to find conditions at an intermediate time at which data 
are desired. 

"There are at least three forms in which the Adams long-werm pre¬ 
dictor-corrector formulas can be written. The only difference is in 
mathematical form; therefore, the accuracy and stability are the same 
for all three. 

"The first form is the conventional one in terms of the successive 
backward differences; the second is in terras of the successive values 
of the function. The third is due to Nordsieck and uses the succes¬ 
sive higher derivatives of the approximating polynomial. Each of the 
three forms has certain advantages and certain disadvantages which 
will be discussed now, 

"The backward difference form is fairly easy to start but inter¬ 
polation is somewhat difficult, and it is virtually impossible to 
change intervals except by use of the starting procedure. The suc¬ 
cessive value form of the method is trivial to start up, but inter- 
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polation involves the Lagrangian interpolation formulas, and changing 
intervals is, again, almost impossible. The inability to change in¬ 
tervals immediately after starting causes, as in the present ITEM pro¬ 
gram, a situation where the starting solution is called 28 times, but 
used only 7 times. 

"The Nordsieck method is fairly difficult to start, but very amen¬ 
able to arbitrary changes of time intervals and to interpolation to 
intermediate points. Five points are all that are needed to start 
after a change in time interval (of about 4:1). Due to its versatility, 
the Nordsieck method of degree 5 (called m ■ 6 by Nordsieck) without 
iteration and without choice of interval is used in the program.” 

Starting Procedure (Runge-Kutta-Gill, RKG. Method) 

r i3~| 

The program uses the Gill modification of Runge-Kutta L J 
for the starting procedure. Runge-Kutta methods are widely 
used and the differences between Gill’s method and others are minor. 
Therefore, only the formulas as given by Gill and not the deriva¬ 
tions are presented. The method was chosen because it introduces some 
simplicity and error reduction over similar methods. 

Given the differential equation 



and the initial values x 0 and y(xo)» y(x 0 +h) is calculated by 



( 122 ) 


where 



yo ■ y(*o) 





( 123 ) 
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The above procedure could be used directly to compute y(x^ + h), 
y(xQ + 2h) , ,,, , etc. However, Gill has introduced a ’’bridging" 
technique between one entry and the next. This modification increases 
the accuracy with little increase in complexity. The new procedure** ^ ^ 

may be summarized as: 

Given x 0 , y 0 » y* Ui) - fCx^ y(x ± )) 

calculate for i - 0 , 1 , 2 , 3 

K i “ a ± (y*(*i) - b ± qi) 
x i+ 1 * x ± + T i h 
yi+i * y(*i+i) - yi + h Ki 



qi+i - q± + 3 

Ki - ci y'(xi) 


1 

a ° " 2 

aj " 1 - 

b 0 “ 2 

f\ *..x 

c 0 -7 

.. — /F 

i 

T ° “ 2 

T i 0 

a 2 " 1 + 

/I b 2 - 1 

C 2 " 1 + f\ 

1 

t 2 - 7 

1 

a 3 - - 

b 3 - 2 

1 

C 3 - 2 

t 3 - 0 

For the 

first entry into the 

procedure, qg " 0 , 

and the proce- 


dure is the same as given in Equations (122) and (123). For subse¬ 
quent entries, qg is set equal to the previous qi* • 

It should be noted that the program must apply the method to the 
three second-order differential equations given by Equation (105) 
(Encke’s method), or by Equation (121) (Cowell’s method). The program 
treats the problem as one of six first order differential equations 
given byj 

yi - fi(t» yi, yz» •••» y6> i “ 1* 2 . 6 

where 


(124) 
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yi(t) - x(t) 

y 4 (t) - x(t) 

y2(t) - y(t) 

ys(t) ■ y(t) 

Y3(t) - z(t) 

y 6 (t) - z(t) 


where (x, y, z) and (x, y, z) are the geocentric position and velo¬ 
city vectors (in Cowell’s method) or the perturbations from the nomi¬ 
nal orbit (in Encke’s method). 

The procedure is programmed as follows: 

1 » i ■ lj 2, •••,6 ^ 

2 • y ■ (yit y 2 # •••» y6> * (R(to)» R(to)) 

3. For k - 1, 2, 3, 4 

A * y - (yi* Y 2 » •••. y6) (R» R) 

B. For i - 1, 2, ..., 6 

1 . K ak (yi - bk ♦ qi) 

2 . yi yi + K • h 

3. qi *■ qi + 3 • K - Ck • yi 

C. R ■*- (yi, y 2 » y3> 

• 

R - (y4. ys* y6> 

t «- t + xk • h 

^ ... 

call DERIV, i.e., cal. R - F(R, R, t) 

• •• 

Hence the program exits with R(to 4* h) , R(to + h) , and R(to + h). 

Also, steps 1 and 2 are only executed upon the first entry into the 

procedure. For second or later entries, the initial value of is 

—► 

the value computed during the last entry for K - 4. Clearly, y 
already contains (R(t), R(t)) and need not be reset. 

Transition from Starting Procedure to Long-Term Integration 

The starting procedure yields the solutions of the six differen¬ 
tial equations and their rates of change at six successive times. It 
is necessary to transform these data into the form required by the 
Nordsieck long-term numerical integration procedure. 
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lor each first order differential equation, the Nordsieck method 
requires the following five higher derivatives evaluated at t - tgJ 


a(t 0 ) 

h y(to) 
2! 


b(t 0 ) 

h 2 yttH) 

(t 0 ) 

3! 


c(t 0 ) 

h 3 ydV) 

(t 0 ) 

A! 


d(t„) 

h 1 * y(V) 

(to) 

5! 


e(t 0 ) 

h 5 y (VI) 

(to) 

6! 



J 

The RKG starting method provides data for y(t) and y(t) at 
the six time intervals up to and including to# i,e,, RKG provides: 


y(to) 


y(t 0 ) 


y(t 0 - 

h) 

rr 

o 

1 

h) 

y(to ~ 

2h) 

y(to - 

2h) 

y(to - 

3h) 

y(to - 

3h) 

y(to ~ 

Ah) 

rr 

o 

1 

Ah) 

y(to - 

5h) 

y(to - 

5h) 


The required values for a(to)# b(to)# c(to)# d(to)# and e(to) 
will be found by using Lagrange f s Interpolation formula to fit a 
power series of degree five to the y(t) data provided by the RKG 
method. The power series will then be successively differentiated to 
obtain the required data. Let 

x - t - t 0 

h 

and let primes denote derivatives with respect to x. Therefore, 


(125) 
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f (t) - h y (t) 
y" (t) - h 2 y(t) 
y'" (t) - h 3 y< lV) (t) 
y""(t) - h 4 y< V >(t) [ 

y""'(t) - h 5 y< VI >(t) 
y"""(t) - h 6 y (VII) (t) 


From Lagrange's Interpolation Formula, 


where 


F 0 (x) 


F 1 (x) 


F 2 (x) 
F 3 (x) 
F 4 (x) 


F 5 (x) 


5 

y(x) - l Fi(x) yi 
i“0 

(x+1)(x+2)(x+3)(x+4)(x+5)^ 

120 

- x(x+2)(x+3)(x+4)(x+5) 

24 

x(x+l) (x+3) (x+4) (x+5) 

12 / 

- x(x+l)(x+2)(x+4)(x+5) 

12 

x(x+l)(x+2)(x+3)(x+5) 

24 

- x(x+l)(x+2)(x+3)(x+4) 

120 J 


and y^ are the values determined by the RKG method. Multiplying 
out the factors in Equation (128) yields 


( 126 ) 


(127) 


(128) 
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Fq (x) 


Fi (x) 


X* + 15. JC. ^ +- m * + 120 

5 j. i/. „4 


120 

. - ^ -H4 x 4 1 71x_ 3 + .13,4 x 2 + 1.2Ox) 

24 

S 4. 107 v 2 


Fz ( x ) . x 5 + .13 x 14 + 59 + 107. x z + 6.0.x 

f 3(x) . - (» 5 «!;»*« 4 | » 3 t_isjsLt «*> 


Ft* (x) 


24 

120 


F 5 (x) 

From Equations (125) and (126) 

*Uo) - ^ 

b(t 0 )-i=5P 

c(to)-^ 


d(t 0 ) 


e(t 0 ) 


4! 

v"" ( 0 ) 

5! 


^. 1*111 
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Successively differentiating Equation (127) with respect to x, 
(using Equation (129)), setting x ■ 0, and substituting into Equa¬ 
tion (130) yields the following in matrix notation; 


2a(t 0 )~ 


"-24 

150 

-400 

600 

-600 

274~ 


*— 

y(to - 

— 

5h) 

3b(t 0 ) 


-50 

305 

-780 

1070 

-770 

225 


y(t 0 - 

4h) 

4c(t 0 ) 

. _L_ . 

-35 

205 

-490 

590 

-355 

85 


y(t 0 - 

3h) 

5 d(t 0 ) 

120 

-10 

55 

-120 

130 

- 70 

15 


*<• 

s 

rt 

o 

1 

2h)| 

_6e(t 0 )_ 


- 1 

5 

- 10 

10 

- 5 

1 


/-s 

rt 

o 

h) 

or 








o 

u 

s-/ 

_i 



( 129 ) 


(130) 
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a(t 0 ) 

b(t 0 ) 

c(t 0 ) 

d(t 0 ) 

^e(t 0 ) 


1 

1440 


-144 

900 

-200 

1220 

-105 

615 

- 24 

132 

- 2 

10 


-2400 

3600 

-3120 

4280 

-1470 

1770 

- 288 

312 

- 20 

20 


-3600 

1644 

-3080 

900 

-1065 

255 

- 168 

36 

- 10 

2 


y(to 

- 5h) 

y(to 

- 4h) 

y(t'o 

- 3h) 

y(to 

- 2 h) 

y(to 

- h) 

y(to) 


(131) 


Equation (131) is used to compute the sets of coefficients [a(to)» 

T 

b(to)> x(to)f d(to)* e(to)] for the six differential equations given 
by Equation (124). 


Long-Term Integration (Nordsieck Method) 

The Nordsieck method is used to continue the solution of Equation (124) 


dy-r 

“ fi(t, ..y6> i " 1* 2 * 

once begun by the RKG method. 

Considering the single equation 


dt 


f(t, y) 


and approximating its solution by a polynomial of degree five 
y p (t 0 +h) - y(t 0 ) + h[y(t 0 ) + rfr y(t 0 ) + y( IIX )(to) 


+ y (Iv ^(t 0 ) + 77 y^ v ^(t 0 ) + 77 y (vx )(t 0 )] 


5! 


6 ! 


(132) 


(133) 
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where 

h is the integration step size, 
to is the value of t at the last integration. 

Substituting from Equation (125) 

y P (t 0 +h) - y(t 0 ) + h[f(t 0 , y) + a(t 0 ) + b(t 0 ) + c(t 0 ) + d(t 0 ) + e(t 0 )] 

Differentiating Equation (133) with respect to t = tg + h and using 
Equation (125) 

f p - f(t 0 , y) + 2a(t 0 ) + 3b(to) + 4c(t 0 ) + 5e(t 0 ) 


Having obtained a predicted value of f(to + h, y(to + b)) a 
value of f (to + h) - f(to + h, y* 5 ) is computed and combined with 
f^ to give the Nordsieck corrector. 


where 


Ki h[f(t 0 + h) - fP] 

„ 19087 

Kl 60480 " °* 315591931 

The values of a(to + h), b(to + h), c(to + h), d(to + h) and 


e(tQ + h) are then computed in terms of their values at t ■ to 

by 


a(to + h) - a(to) + 3b(t 0 ) + 6c(t 0 ) + 10d(to) + 15e(to) 

+ K 2 [f(t 0 + h) - f P ] 

b(to + h) ■ b(to) + 4c(to) + 10d(to) + 20e (to) 

+ K 3 [f(t 0 + h) - f p ] 

c(to + h) - c(t 0 ) + 5d(t 0 ) + 15e(t 0 ) + K 4 [f(t 0 + h) - f p ] 

d(t 0 + h) - d(t 0 ) + 6e(t 0 ) + K 5 [f(t 0 + h) - fP] 

e(t 0 + h) - e(t 0 ) + K 6 [f(t 0 + h) - fP] 

where 

k 2 - “ 1.141666667 


K 3 - | - 0.625 

K 4 - 77* - 0.1770833333 
yo 


( 134 ) 


(135) 


(136) 
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0.025 


k 5 

k 6 


-L_ 

40 

- 0.00138888889 


The Nordsieck method may now be summarized as follows. For each 

p p 

entry, y p is computed by Equation (134), f by Equation (135), and 
f(tQ + h) - f(tQ + h, yP) by the functional relationship implied by 
Equation (132), Then y(t Q + h) is found by 

y(t 0 + h) ■ y^ + Ki h(f(tQ + h) - f p )• 


Finally the coefficients a(t) through e(t) are updated by Equa¬ 
tion (136). 


The integration interval, h, is readily changed; the change 
being accomplished by using new values of a(t), b(t), c(t), d(t), 
and e(t). These new values are obtained from the following equations: 

hn 
h 0 


B 


a n (t) - Ba 0 (t) 
b n (t) - B 2 b 0 (t) 
c n (t) - B 3 c q ( t) 
d n (t) - B 4 d Q (t) 
e n (t) - B 5 e 0 (t) 


where the subscripts n and 0 stand for new and old, respectively. 
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COORDINATE SYSTEMS AND TRANSFORMATIONS 

The equations of motion of a vehicle are given by Equation (2). 

It is necessary to express the vectors used in this equation in an 
inertial coordinate frame, i.e., a coordinate system in which inertial 
forces due to the system^ acceleration are absent or at least negli¬ 
gible. A coordinate system rigidly attached to the Earth is inappro¬ 
priate, since such a system experiences noticeable accelerations due 
to the daily rotation of the Earth. 

The inertial frame chosen for this program is a so-called "mean 
equinox of base date” system which is determined by the vernal equi¬ 
nox of 0.0 hours 1 January of the year subsequent to the epoch time 
(initial input time) used in a given run. This particular base date 
system has been chosen as a basis for calculation because the plane¬ 
tary, lunar, and solar coordinates are written on tapes in that coor¬ 
dinate system. Rather than transform the tape information, the vehi¬ 
cle initial conditions and the oblateness accelerations are transformed 
into the base date system. 

The definition of the base date system and various aurxliary 
Earth referenced systems is given below along with the transformations 
between them. 

Terminology 

First, consider the celestial sphere (Figure 14) which is a 
sphere of infinite radius with the Earth at its center upon which the 
positions of stars, planets, the Sun, and the Moon are projected. It 
is sometimes convenient to pick the center of the sphere to be at an 
observer or the center of the Sun, but for the purpose of our discus¬ 
sion its center will always be taken to be at the geocenter. The 
stars may be taken to be fixed on the sphere (their small proper 
motions being neglected) while the Sun, planets, and the Moon appear 
to move and describe certain paths along its inner surface. The ap¬ 
parent path of the Sun during the course of a year on the celestial 
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sphere is a great circle called the ecliptic . The north celestial 
Pole is the point where the extension of the Earth’s north pole inter¬ 
sects the celestial sphere and the south celestial pole is defined 
similarly. The great circle where the extension of the Earth’s equa¬ 
tor intersects the celestial sphere is called the celestial equator . 

The point on the celestial equator where the apparent path of the sun 
crosses into the northern half of the celestial sphere at the begin¬ 
ning of spring is called the vernal equinox . This is one of two 
points where the ecliptic and equator intersect. In the following 
discussion the word celestial will usually be dropped and we will 
just speak of the north pole or the equator. 

The equator and ecliptic are constantly in motion due to the ef¬ 
fects of nutation and precession. It should be pointed out that this 
is astronomical precession which is distinct from the force free pre¬ 
cession due to the flattening at the poles. The former is caused by 
the net torque on the equatorial "bulges** due to the gravitational 
attraction of the sun and moon. The torque is quite small so that the 
precession is extremely slow - the period being 26,000 years compared 
to the rotational period of one day. The total applied torque is not 
constant in time, because the torques of the Sun and Moon have slightly 
different directions to the ecliptic and vary as the three bodies move 
around each other. As a result, there are irregularities in the pre¬ 
cession, commonly designated as astronomical nutation . The astronomi¬ 
cal nutation is not to be confused with the "true nutation** which is 
present even in the force-free precession of the Earth’s rotation 
axis. In the subsequent discussion the astronomical precession and 
nutation will simply be referred to as precession and nutation. 

Precession may be separated into luni-solar precession, and 
planetary precession. The luni-solar precession causes the north 
pole to move around the ecliptic pole . The ecliptic pole is the 
point of intersection of the extension of the normal to the ecliptic 
plane at the geocenter, and the celestial sphere. Nutation is an 
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NORTH CELESTIAL POLE 



ECLIPTIC 


CELESTIAL EQUATOR 


Figure 14. Celestial Sphere 
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Figure 15. Precession 
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irregular motion of the pole with a period of about 18.6 years. These 
two effects combine to model the motion of the true equator . The ac¬ 
tual equator at any time is called the true equator of date. The word 
date here refers to any arbitrary time. The mean equator of data i« 
defined by the location of the equator when the effects of nutation 
are removed (i.e., where the true equator of date would be if there 
were no nutation). 

The gravitational action of the planets on the Earth causes the 
plane of the Earth’s orbit to slowly precess. This motion appears 
from the earth as a slow precession of the ecliptic, and is called 
planetary precession. This effect causes the obliquity (the angle 
between the equator and the ecliptic) to decrease at the rate of 
about 47" a century. Therefore, to be precise we refer to the eclip¬ 
tic of a specified date. 

The effect of precession is illustrated in Figure 15. The mean 
pole, mean equator, and ecliptic at time to are labeled base date 
and at a later time, t, are labeled date. Luni-solar precession 
causes the mean equator to move between base date and date, while 
planetary precession moves the ecliptic. Figure 16 illustrates how 
nutation causes the true equator of date to differ from the mean 
equator of date. 

In this discussion we will use several coordinate systems which 
will be defined in terms of the tyue equator and ecliptic and mean 
equator and ecliptic at a given time. The true vernal equinox (or 
true equinox) is defined by the point of intersection of the true 
equator and ecliptic. The mean vernal equinox (or mean equinox) is 
defined by the point of intersection of the mean equator and the 
ecliptic. The true equinox and the mean equinox are points on the 
celestial sphere that experience a slow motion and hence we must 
specify a time to indicate their exact positions. 
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Figure 16. Nutation 
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Coordinate Systems 

The five coordinate systems used by the program for generating 
trajectories are: 

(1) The mean equinox of base date system * This coordinate 
system is tne inertial system employed by the program 
and employs unit basis vectors x, y, and z defined 
as follows: 

x is a unit vector directed towards the mean vernal 
equinox of base date* 

z is a unit vector normal to the mean equatorial 
plane of base date, positive in the northern 
hemisphere. 

y is a unit vector orthogonal to x and z. 

(2) The mean equinox of date system* This system employs 

unit basis vectors x^ f yy, and zm defined as 
follows: 

x^f is a unit vector directed toward the mean 
vernal equinox of date. 

ZH is a unit vector normal to the mean equatorial 
plane of date. 

y^l is a unit vector orthogonal to xm and z^. 

(3) The true equinox of date system . This system 

employs unit basis vectors x-p» yii and zx 
defined as follows: 

xx is a unit vector directed toward the true 
vernal equinox of date. 

Zj is a unit vector normal to the true equatorial 
plane of date. 

yx is a unit vector orthogonal to xx and zj» 

(4) The Greenwich coordinate system * This system 

employs unit basis vectors xg» yG» and Z G 
defined as follows: 

xq is a unit vector directed toward the inter¬ 
section of the Greenwich meridian with the 
true equator of date. 

zq is a unit vector normal to the true equator 
of date* 
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yG is a unit vector orthogonal to xq and zq . 

(5) The topocentric coordinate system . This system employs 

three orthogonal basis vectors x g , y« and z g , 
directed toward the east, north, and local vertical, 
respectively, at a point on the Earth's surface 
(considered as an ellipsoid)* It is only used for 
observation calculations and its details are dis¬ 
cussed under OBSERVATIONS. 

The differences between the above coordinate systems are due to 
dynamical effects of the Earth's motion or to the geometry of the 
Earth** It is possible to express a 3 * 1 vector of position or velo¬ 
city written in the basis of one coordinate system in terms of another 
coordinate basis by applying 3 * 3 orthogonal matrix transformations. 

—f 

Table III below lists the transformations used by the program. R is a 
3 x 1 vector whose subscript indicates the basis in which the vector 
is expressed* 

Table III. 

Transformation Matrices. 


MATRIX 

EFFECT 

FROM 


TO 


EQUATION 

[G] 

Earth Geometry 

*s 

ys 

*s 

X G 

yG 

ZG 

R G - [G] R s 

ty] 

Right Ascension of 
Greenwich from True 







“ tY] Pg 

Equinox of Date 

*G 

yg 

ZG 

x T 

y T 

zx 

[N] 

(Daily Rotation) 

Nutation 

Xj 

H 
1 >> 

Z T 

xm 

yM 

ZM 

“ IN] R T 

[P] 

Precession 

x M 

yM 

z M 

X 

y 

z 

R - [P] R M 

The derivations of [y]. 

[N] 

• 

and 

[P] 

are 

given below. 


should be recognized that these three matrices are functions of time. 
Sidereal Time 

When computing the inertial coordinates of a point or station 
on the surface of the rotating Earth at a particular time, a relation¬ 
ship must be found between the Greenwich meridian (the great circle 
that passes through the true north pole and Greenwich) and the vernal 
equinox at that instant. The angle measured in the equatorial plane 


77 












from the Greenwich meridian to the true vernal equinox is called the 

apparent Greenwich sidereal time . 

The right ascension of the mean equinox measured from the true 

equinox is called the equation of the equinoxes or nutation in right 

ascension and is usually denoted by 6a, The relationship between 

the apparent sidereal time and the mean sidereal time is: 

Mean sidereal time - Apparent sidereal time - 6a, 

6a is always measured along the true equator from the true equinox 

eastward. The apparent and mean Greenwich sidereal times are tabu- 

[41 . 

lated in the American Ephemeris and Nautical Almanac J for 0 n of 
each day of the year along with the equation of the equinoxes, 6a. 
Each entry is given in hours, minutes and seconds of time to 0 s .001, 
where 

l h - 15°, l m - 15', and I s - 15". 

Apparent Greenwich sidereal time, Oq* since it is measured 
from the true equinox, is affected by nutation and hence changes at 
an irregular rate (see Figure 17), Since mean sidereal time, 0 q, 
is measured from the mean equinox, which is affected only by preces¬ 
sion, it increases at almost a constant rate with only a small accel¬ 
eration and may be calculated at 0* 1 of any day by the formula 
0G - 6 h 38® 45?836 + 8,640,184?542T + 0?0929T 2 

where T is the number of Julian centuries (36,525 days) from noon 
0 January, 1900 (Julian day 2,415,020.). 

An equivalent formula in degrees is 

0 G (d) - 100!0755426 + 0!985647346 d + 2S9015 * 10" 13 d 2 

where d is the integer number of days past 0^ 1 January, 1950. 

Equations (137) and (138) may be reformulated as a function of the 
time past any conveniently chosen epoch resulting in an equation of 
the same form with different coefficients. 

Evaluating the above formula for 0 G (d + 1) and 0 G (d), we 


(137) 


(138) 
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can find the angle through which the Greenwich meridian rotates in a 
day with respect to the moving mean equinox. 

0G(d + 1) - 0G(d) + 360? - 360?985647346 + 5? * 10 -13 d. 

in which the additive term 2?9015 * 10 _13 has been neglected) Di- 
viding Equation (139) by the number of seconds in a day (86,400) and 
neglecting the final term whose contribution is insignificant over a 
period of only a few years ( < 25 yrs.), we have w s , the rotational 
rate of the Earth with respect to the moving mean equinox, 

u B - 0?0041780 74622 sec. -1 

Therefore, to find the mean Greenwich sidereal time at time 
d + t we use 

©G(d + t) - 0G(d) + u s t 

where d represents the integer number of days past 0* 1 1 January, 

1950, and x the number of seconds that have elapsed since 0^ of 
the day. 

The apparent Greenwich sidereal time 0 q can be found at time 

t by 

6 € <t) - 0G(t) + 6a. 

Now, w s will be referred to as the sidereal rotational rate 
of the Earth. The word sidereal is used since the rate is with res¬ 
pect to the mean equinox, which is the reference point used for 
sidereal time. 

Transformations 

The effects of nutation and precession are fairly small and 
over a period of a few days the errors introduced by neglecting these 
effects is on the order of ten meters at most. Therefore, for short 
trajectory determination runs the program has the option of not in¬ 
cluding them by setting the [P] and [N] matrices equal to the 
identity matrix, [I]. The [y] matrix which includes the effects 
of daily rotation must be included, however. 


( 139 ) 


(140) 


(141) 


(142) 
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Figure 17. Sidereal Time 
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Referring to Figure 17 the transformation from the Greenwich co¬ 
ordinate system to the true equinox of date system consists of a ro¬ 
tation along the true equator through the angle 0c(t) where 0G(t) 
is given by Equation (142). Thus the [y] matrix is given by: 


[y] 


cos 0c(t) - sin 0 g(O 0 

sin 0 g( t) cos 0^(0 0 

0 0 1 


(143) 


If the effects of nutation and precession are to be neglected then 
0G(t) replaces 0^(0 in Equation (143), and the reference system 
of the computation may be considered as the system defined by the 
mean equinox of date (i.e., the effects of nutation and precession 
are neglected), 

Next, an expression for the nutation matrix [N] is found. 

In Figure 18, y is the true equinox of date and y is the 
mean equinox of date. The right ascension of the mean vernal equinox 
referred to the true equinox is 6a, the equation of the equinoxes. 
Notice that 6a is less than zero as it is drawn in Figure 18, since 
right ascension is measured positively toward the east. The nutation 
in longitude, A^, is the longitude of the mean equinox measured 
from the true equinox. Celestial longitudes are measured along the 
ecliptic, and since the positive direction is east, A^ is also neg¬ 
ative in Figure 18. The mean obliquity, e, is the inclination of 
the ecliptic to the mean equator, and e, the true obliquity, is the 
inclination of the ecliptic to the true equator. The nutation in 
obliquity is 6e, where 

6e - e - e. 


From inspection of Figure 18, and by using Napier’s rules for 
right spherical triangles, 

sin(90° - e) - tan 6a tan(90° - A^) 
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or 


tan 6a - tan Aij; cos e 


Using the series expansion for the tangent function we have 


6a + 

Now, 

fore 


(6a)3 2 C6o) 5 

3 15 * 

| <5a |, and | 6<|> | 

, Che approximation 


are less than 1* * 10~ 4 radians. There- 


6a - A\J/ cos c 


should never produce an error greater than 1 * 10~ 12 radians or 
2 * 10 7 seconds of arc. 

The nutation matrix [N], is now derived which transforms co¬ 
ordinates referred to the true equinox of date into coordinates re¬ 
ferred to the mean equinox of date. This transformation consists of 
three rotation matrices. Using Figure 18, we note that the first, 
[A], is a rotation about the x-axis through the angle e. The sec¬ 
ond, [B], is a rotation about the z’-axis through AiJ/, and the 
third, [C], is a rotation about the x"-axis through - e. 
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[N] - [C][B][A] 


[C] - 


10 0 
0 cos e - sin e 

0 sin e cos e 


cos 

- cos e sin 

- sin e sin A\Jj 


sin A^ cos e 

cos e cos A^ cos e 
+ sin e sin e 

sin e cos A^ cos e 
- cos e sin e 


sin A\[» sin e 

cos e cos A^ sin e 
- sin e cos e 

sin c cos A\J) sin e 
cos e cos e 


(144) 


Since || < 10" 14 and |e - e| < 10“^ the nutation matrix is 
often approximated by the following matrix; 


[N 1 ] 


1 

- A^ cos e 

- Aip sin e 


Av cos e 
1 

- 6c 


A^ sin c 
6e 

1 


The error in the elements in [N 1 ] should be less than one unit 
in the eighth significant figure^ 

Reference [ 14] gives formulas for A^ and 6e which depend on 
the mean longitude of the Sun, the mean longitude of the perigee of 
the Sun, the mean longitude of the Moon, the mean longitude of the 
ascending node of the Moon, and the mean longitude of the lunar peri¬ 
gee, There are 69 terms in the expression for A^ and AO terms in 
6e, Truncated expressions for and 6e are given in Reference [l]. 

Finally, an expression for [P], the precession matrix is de¬ 
rived below. Figure 19 illustrates the effect of precession over the 
period of time t - to* In this figure, 90° - * s t ^ ie right ascen¬ 
sion of the ascending node of the mean equator at time t on the 
mean equator at to measured from the mean equinox of to; 90° + z 
is the right ascension of the node measured from the mean equinox of 
t; and 6 is the inclination of the mean equator at time t to the 
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Figure 19. Formation of Precession Matrix 
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mean equator of to# 

The precession matrix, [P] , is now derived which transforms 
coordinates referred to y(t) to coordinates referred to y(to)# 

[P] is a result of three rotation matrices. From Figure 19 we 
see that the first, [A 1 ], is a rotation of (90° + z) about the z- 
axis. The second, [B 1 ], is a rotation about x through -0 and 


the 

third, [C’]. 

is 

rotation about 

z 

through (Co 

- 90°). 


[P] - [A’] 

[B 1 ] 

[c'J. 

where 







COS 

/-s 

VO 

o 

0 

+ z) 

sin 

(90° + z) 

0 


r» 

- sin z 

cos z 

[A'] 

- 

- sin 

(90° 

+ z) 

cos 

(90° + z) 

0 

- 

- cos z 

- sin z 



0 


0 


1 


0 

0 




1 


0 

0 


1 

0 

** 

0 


[ B f ] - 

0 

cos 

(-0) 

sin (-0) 

- 

0 

cos 0 

- sin 0 




0 

- sin 

(-0) 

cos (-0) 


0 

sin 0 

cos 0 


[C'J 


cos (co - 90°) 
- sin (Co - 90°) 
0 


sin (; 0 - 90°) 0 

cos (Co ~ 90°) 0 

0 1 


sin Co ~ cos Co 0 

cos Co sin Co 0 
0 0 1 


Hence, 



sin z 

sin Co 



COS z 

sin 

40 

sin 

0 

cos 

40 


+ cos 

z cos 0 

cos 

40 

+ sin 

z cos 0 cos Co 





- 

- sin z 

cos Co 



cos z 

cos 

40 

- sin 

0 

sin 

40 


- cos 

z cos 0 

sin 

40 

- sin z 

cos 

0 sin Co 






- cos z 

sin 0 



- sin z 

sin 

0 

cos 

0 




Reference [ 14] gives the following numerical expression for Co» a » 
and 0 for the case when coordinates are precessed to a later epoch. 
Initial epoch, to ■ 1900.0 + Tq 
F inal epoch, t - 1900.0 + tq + t 

C 0 * (2304."250 + l."396 t 0 ) t + 0."302 t 2 + 0,"018 t 3 
z - Co + 0."791 t 2 

0 - (2004."682 - 0,"853 t 0 ) t - 0."426 t 2 - 0."042 t 3 


(145) 
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where t 0 and t are in units of tropical centuries. 

When the initial epoch is 1900.0 + Tq + t and the final epoch is 
1900.0 + tq (i.e., we are precessing a vector to an earlier epoch) 
one whould replace 
Co by - 2 
z by - t 0 
0 by - 6 

in the above matrices [A 1 ] , [B 1 ], [C 1 ], and [P]* 

Combining the results of Equations (143), (144), and (145) with 
the relations in Table III, all necessary coordinate transformations may 
be performed by the program. 

OBSERVATIONS 

The SPACE-A program has a provision for a total of 25 observation 
types of which a number of ground-based observations are specified and 
described here. The specified observations include: 


azimuth, A 

elevation, E 

round-trip range, 2p 

. 

one-way range rate, p 

hour angle, H 

declination, D 

il-direction cosine, i 

m-direction cosine, m 

X-angle, X 

Y-angle, Y 

range equivalent time, At 

range-rate equivalent time, At* 


The position and velocity of a vehicle is available at all times 
from the trajectory generation portion of the program. The station 
position and velocity is computed from input data concerning its geo¬ 
detic coordinates. Therefore, one may consider the relative geometry 
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of the vehicle and station to obtain the observations. 


The program has the option to include the effects of refraction 
upon ground-based observations. A simple model of the Earth’s atmos¬ 
pheric effects upon observations is employed in which the index of re¬ 
fraction varies only with altitude. Neglecting the effects of refrac¬ 
tion causes the predicted elevation angle to be less than the observed 
elevation. Errors are also made in the computation of range and range 
rate related observations. 

The method used to include the effects of refraction is to com¬ 
pute range and range rate from the relative geometry of the station 
and the actual vehicle position and velocity and then to. apply correc¬ 
tions based on the refraction model in order to obtain the range and 
range rate related observations. In the case of angle related mea¬ 
surements a slightly different approach is taken. A vector pointing 
from the station to the computed vehicle position is constructed. 

Then the refraction model is used to obtain a new vector pointing from 
the station to the observed vahicle position. All angle observations 
are then computed directly or indirectly using this new vector. 

Coordinates Used in Observation Calculations 

The input coordinates of the station (the station is assumed to 
be on an ellipsoidal Earth) are initially given as geodetic longitude, 
(X G ) t geodetic latitude (<J>q) , and height. Given this information 
the position and velocity of the station in the true system of date 
coordinate system is given by: 

x T - [c + h] cos 4 >q cos (Xq + Y) 

y x ■ [c + h] cos 4 >g sin (X G + y) 

z T - [c(1-e 2 ) + h] sin <J> G V 

x x - - w e y T 

y T - u> e x x 

z X ■ 0 


( 146 ) 
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C “ 



*e 

e 2 sin 2 $g 


(U7) 


e 2 - 2f - f 2 

where 


(148) 


X T* yT* zt are the true system of date coordinates of 
the station, 

Ag, 4>g» h are the geodetic coordinates of the station, 

y is the right ascension of the Greenwich 
meridian in the true system of date, 

R e is the equatorial radius of the Earth, 

, i 

f is the flattening constant of the Earth 


and -- --*- - '-298.30- 

Using the precession matrix, [P] and the nutation matrix, [N], we 

can express the station coordinates in the inertial (base date) coor¬ 
dinate system. From the position and velocity vectors of the station 
and the vehicle one can obtain the vectors of position and velocity of 
the vehicle with respect to the station expressed in the inertial sys¬ 
tem. Thus 

P - P s - R 

i - L - i 


(149) 

(150) 


where p, p are vectors of the vehicle with respect to the 
station, 

Pat P® are the vectors of the station, and 
R, ^ are the vectors of the vehicle. 

Note that all of these vectors are expressed in the inertial coordi¬ 
nate system. 

For observation calculations it is convenient to define a station 
oriented coordinate system which has a basis consisting of three unit 
orthogonal vectors pointed toward the east, north and local vertical. 
Using station oriented coordinates the unit vectors are: 
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Xg - [10 

o] T , 

east vector 

ys ■ [° i 

o] T , 

north vector 

zs - [o o 

i] T , 

up vector 


Expressed in terms of the inertial base date system these vectors are: 



e - [T] x s 

■> , , 4 - 

east vector 

(151) 


n - [T] y s 

north vector 

(152) 


u - [T] z S 

up vector 

(153) 

and 

[T] - [P][N][ Y ][G][o] 

(154) 


where all vectors are 3 x 1 vectors. 


[P],[N], and 


[a] is a 3 x 3 rotation matrix to account for 
misalignments between the station oriented 
coordinates and the true topocentric east, 
north, and vertical coordinates, 

[G] is a 3 x 3 matrix transforming coordinates 
expressed in the topocentric system into 
coordinates of the Greenwich system, 

[Y] are the 3 x 3 matrices of precession, nuta¬ 
tion, and right ascension of the Greenwich 
meridian which transform the coordinates into 
the inertial base date system. 


[G] 


sin Xq - sin $q cos Xq 

cob Xq - sin $G sin Xq 

0 cos <Pq 


cos 4 >g cos ^G 
cos <|>g sin Xq 
sin <J>g 


055) 


where Xq and <Pq are, respectively, geodetic longitude and geo¬ 
detic latitude of the station. 

The final result of the above transformations is that the station 
oriented basis vectors (e, n, and u) are expressed in the inertial 
base date coordinate system. 


Computation of Observation Types 

In the following discussion all computations are made using vec¬ 
tors expressed in the inertial (base date) coordinate system. 

Azimuth and elevation measurements are illustrated in Figure 20, 
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Azimuth, A, is measured easterly from station north from 0 £ A £ 2n. 
Elevation, E, is measured from the stations horizontal plane upward 

TT TT “► 

with a range " ” i E i “< If P is the station to vehicle vector 
the azimuth and elevation are given by: 


A - 



E • tan 


-1 


/ 


0 * U 


(p • e) 2 + (p • n)' 


If refraction is to be included in the computations, the vector, p, 
is replaced by the vector p', where p' is the vector from the 
station to the observed vehicle position as described below. 


(156) 


(157) 



Figure 20. Azimuth and Elevation 


The round-trip range, 2p, is twice the distance from the station 
to the vehicle and its value is given by: 

2P - 2|p| + 2Ap (158) 

If the effects of refraction are not included in the calculation 
2Ap - 0. If the refraction effects are included Ap is a correction 
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terra the value of which is given by Equation (202). 

The one-way range-yate , p, is equal to the magnitude of the 
time derivative of the vector from the station to the vehicle and its 
value is given by: * 


P 



4* Ap 


( 150 ) 


If the effects of refraction are not included Ap - 0. If the effects 
of refraction are included Ap is a correction terra the value of 
which is given by Equation (206)• 

The hour angle , H # is the angle between the station meridian 
and the projection on the true equator of the station to vehicle vec¬ 
tor measured in the true equatorial plane. It is measured westward 
from - tt ^ H <_ tt. The declination . D, is the angle made with the 
true equatorial plane by the station to vehicle vector. Declination 
is measured positive in the northern hemisphere and has a range 
- — <_ D <. Hour angle and declination are illustrated in Figure 21 
which shows the station at the center of a celestial sphere and the 
projected true equator on the celestial sphere. The hour angle and 
declination may be derived using spherical trigonometric relations in 
terms of geodetic latitude, azimuth and elevation (<$, A, E), The 
values are given by: 


H 


tan" 1 


sin A 

cos A sin 4 >q - cos tan E 


(160) 


D - sin" 1 [sin 4 >q sin E + cos 4 >q cos E cos A] 


(161) 


The values of A and E used above are given by Equations (156) and 
(157). 

The ^-direction cosine . £, and m-dire,ction cosine . m, are 

shown in Figure 22. The ^-direction cosine is the cosine of the angle 

-¥ 

between tne station to vehicle vector p, and the station east vector, 
—► 

e. The m-direction cosine is the cosine of the angle between the sta- 
tion to vehicle vector p, and the station's north vector, n. The 
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values of £ and m are given by: 



( 162 ) 


m - 



(163) 


If refraction effects are included p f replaces p in the above 
equations. 

The X-angle and Y-angle measurements are illustrated in Figure 23. 
The Y-angle is the angle between the station to vehicle vector, p, 
and the perpendicular projection of this vector on the station’s east- 
vertical plane. It is positively measured toward the north, negatively 
toward the south with limits - ~ <. Y The X-angle is measured 

between the vertical vector, u, and the perpendicular projection of 
the station to vehicle vector onto the station’s east-vertical plane. 

It is measured from the positive vertical eastward with limit 
” * .1 X The value of X and Y are given by: 


X - 



Y - tan 


-1 


/ 


P » n 


(p • e ) 2 + (p • u ) 4 


If refraction effects are included p’ replaces p in the above 
equations. 

The range time equivalent . At, and the range-rate time equiv¬ 
alents At’, are included since the raw data from typical tracking 
systems are the time between a transmitted and received signal for 
range, and the time to count a given number of Doppler cycles for 
range-rate. In most systems, these quantities are first converted to 
range and range-rate units. However, it may be found useful in some 
cases to use the raw data equivalents. The values of At and At' 
are given by: 

At ■ 

c 



(164) 


(165) 
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Zenith 



Figure 21. Geometry Illustrating Hour Angle and Declination 



Figure 22. i!-cosine and m-cosine Figure 23. X-Angle and Y-Angle 
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At* - 


(167) 



c 


where 

c is the velocity of light, 

|Pl| and |p 2 | are the ranges at the beginning and end 
of the measurements, respectively. 

Specific provision is made in the program for computing doppler 
shift frequency measurements. The equations for these observations 
are not included, since the method of computation can change depending 
on the specific hardware in use for a particular mission. 

Effects of Refraction 

Because of variations in the refractive index of the atmosphere 
the propagation path of a tracking beam is subject to refractive 
bending. 

Two models for the variations of refractive index are included 
in the program: one is a model for the troposphere, the other is a 
model for ionosphere, A numerical method due to Weisbrod and 
Anderson^is employed in which the effects of refractive bending 
are determined by numerically integrating over the total propagation 
path, the index of refraction at each point being determined by the 
assumed model. The vector from the station to the current calculated 
position of the vehicle, p, is available from the program. Using 
numerical methods a correction to the elevation angle is found and a 
new vector, p 1 , from the station toward the observed position of 
the vehicle is determined. 

In addition to angle corrections due to refractive bsnding, the 
effects on range measurements due to signal retardation, and the ef¬ 
fects on range-rate measurements due to refractive bending are also 
found• 

Index of Refraction Models 

In order to simplify computational problems, atmospheric models 
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representative of average conditions are employed in which the fol¬ 
lowing assumptions are made: 

(1) The gradient of the index of refraction varies only 
with altitude, 

(2) The index of refraction profile is approximated by a 
number of linear segments in which the length of each 
segment is very small compared to the Earth’s radius, 

( 3 ) The troposphere extends to 40 kilometers, 

(4) The region between the end of the troposphere and the 
beginning of the ionosphere is assumed to have zero 
refractivity (and therefore no bending or signal re¬ 
tardation occurs), 

( 5 ) The ionosphere lies between a height ho (input data) 
and 2,000 kilometers, 

(6) The index of refraction is zero beyond 2,000 kilometers. 

In the tropospheric model, refractivity (N) is assumed to decay 

exponentially, with the ground index of refraction and the scale 
height as parameters. The equation for the tropospheric model is as 
follows: 

N - N 0 e _h/H - (n - 1)10 6 

where 

Nq is the refractivity at sea level, an input quantity 
for each station (usually 313), 

h is the height above the Earth, 

H is the scale height, an input quantity for each 
station (usually 7 kilometers), 

n is the index of refraction. 

For the tropospheric model, the refractive errors are considered 
to be independent of signal frequency since the index of refraction 
is virtually independent of frequency up to 30,000 MHz, 

In the ionospheric model, the index of refraction is dependent 
upon more parameters than those considered for the tropospheric model. 
The ionosphere consists of several belts of charged particles. The 
F layer is very much larger than any other layer, and contains a 
greater number of charged particles than the other layers. The F 


( 168 ) 
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layer is the one closest to the Earth’s surface. It is subdivided 
into the FI and F2 layers. In the ionospheric model, the index 
cf refraction is primarily dependent upon the height, ho, of the 
base, of the ionosphere’s F2 layer, the maximum electron density of 
the F2 layer, and the height of the maximum electron density of the 
F2 layer. 

Both index of refraction and the height, ho* are dependent upon 
diurnal, solar activity, seasonal, and geographical variations as well 
as other miscellaneous sporadic variations. Unlike the tropospheric 
model, the refractive errors in the ionospheric model are frequency 
dependent. 

In constructing the model, the range of the signal frequency has 
been limited to frequencies above 100 MHz, since this range of spec¬ 
trum both represents the situation of greatest interest and enables 
equation simplification. 

The relationship between the index of refraction (n) the angular 
frequency of the incident signal (w), and the electron density in the 



(169) 


where 

P e is the electron density per cubic meter, 

e is the electron charge (1,6 * 10~ 19 Couloumbs), 

m is the electron mass (9,08 * 10~ 31 kilograms) 

is the permittivity of free space (8,854 x 10~ 12 
farads/meter)• 

Using the first two terms of the binomial expansion as an approx¬ 
imation, the equation for the index of refraction reduces to: 


n - 1 - 40,3 10~ 1? 

Note that f has the dimensions of MHz, This 

L TT 


( 170 ) 


where f - 


equation is true for frequencies above the critical frequency, f c , 
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which is defined as: 


f c - 8.97 po ^ * 10~ 6 (MHz) 

where pg is the maximum electron density per cubic meter. 

From the definition of N in Equation (168), Equation (170) can 
be written as 

p e -c 
N - - 4.03 ~ * 10 5 

The modal selected for electron density versus height consists 
of a parabolic variation below the height of maximum electron density 
matched to a hyperbolic secant profile above the maximum. The rela¬ 
tionships are as follows: 

Pe - Po tl - U - o) 2 ] for 0 <, o <, ll 

p e ■ sach [•“ (o - 1)J for o > 1 J 

where 

h - ho 
° " hm - h 0 

h is the height above the Earth, 

hg ia the height of the base of the F2 layer (input 
quantity), 

h m is the height of the maximum electron density in 
the F2 layer (input quantity). 

Figure 24 ia a plot of the ionosphere model normalized with res¬ 
pect to o and 1/2 (p e /pg). The h, t^, pq parameters refer to 
the ionosphere's F layer. Using this model, the refractive effects 
of the D and E layers are not singled out, because they are quite 
small in comparison with those due to the F layer and are approxi¬ 
mately accounted for by allowing the electron density at the bottom 
edge of the F layer to be zero. 


(171) 


(172) 


(173) 
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Figure 24. Normalized 3-Parameter Model of Atmosphere 













Computation of Ray Bending 

The major effect of ray bending due to refraction is that the 
observed elevation angle, E, of a vehicle in view of a station is 
greater that the actual or computed value of elevation, E c . This 
effect is clearly seen in Figure 25* The difference between the com¬ 
puted and observed elevation is the angle 6. 

In Figure 25 three vectors of importance are shown: 

—► 

p, the vector from the station to the vehicle, 

P 1 , a vector from the station toward the observed 
vehicle position (i,e,, tangent to the ray 
path at the station), and 

Ap, a vector constructed perpendicular to p. 

The value of azimuth. A, and computed elevation, E c , are ob¬ 
tained from Equations (156) and (157) using the value of p, since 
these values do not incorporate refraction effects. From the geome- 
try of Figure 25, one can find the components of Ap in the east, 
north, and vertical station-oriented coordinate system. Thus: 


( 174 ) 


( 175 ) 


( 176 ) 


Ap ■ |l?| tan 6 [T] 


sin E c sin A 
sin E c cos A 
- cos E c 


where |p| is the magnitude of the vector p, and [T] is the 3 * 3 
matrix of Equation (154) used to transform the coordinates of Ad 
from the station oriented coordinate system into the inertial system. 
From the figure it is also obvious that 

p + Ap * p 


or 


■i . 


p - 1^1 tan 6[T] 


sin E c sin A 
sin E c cos A 
- cos E„ 


In the program the vector p 1 is actually normalized to a unit vec¬ 
tor. The length of p 1 is not important, however, since it is used 
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Figure 25. Geometry of Ray Path Used to Obtain 
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only to calculate observed angle-related measurements such a9 azimuth* 
elevation, ^-cosine, m-cosine, X-angle, Y-angle, and indirectly, hour 
angle and declination as described above. 

In Equation (176) all quantities are easily computed with the 
exception of 6, the "error" in the elevation angle. The method of 
computing 6 used by the program is that due to Weisbrod and Anderson*" . 

Only an outline of their analysis is presented here. 

First, consider the ray of Figure 26 entering an infinitesimal 
layer of thickness dr at an angle 6. At a boundary between two 
infinitesimal layers Snell’s law holds. Thus: 

n cos 3 * constant (177) 

where n is the index of refraction. Differentiating Equation (177) 
with respect to path length, s, yields Equation (178) below. Since 
n varies only with altitude or distance along an earth radius, r, 
one obtains: 



4^ cos 6 - n sin 6*0 

as as 


(178) 


dn dr . d6 . _ n 

— i- COS 6 - n 37 sin 6 - 0 


(1/9) 

but 

- sin B, 
ds 9 


(180) 

and from the definition of curvature one gets: 




m 1 

ds K 


(181) 

where < is the 

radius of curvature of the ray. 

Thus, after substi- 


tution, Equation 

(179) becomes: 




1 1 dn 0 

— - — cos 6 . 

k n dr 


(182) 

From Figure 

26, an infinitesimal length ds 

of the ray path is 


given by: 

dr 

sin 6 


(183) 
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thus, 


dY 


^ sin 0 * 


Combining Equations (184) and (182) raaults in an expression for 


dyi 


dY • * ^ cot 8 dr 

Also, from Figure 26, it is aeon that th« dy'a of all elementary 
layers are directly additive. Therefore, the angle Yjk due to the 
total bending between layers bounded by the heights rj and r^ is 
simply.the integral of Equation (185): 


(184) 


(185) 


Yjk 


rit 

7 cot 6 dr 
n or 


r j 

If a ray departs from a given layer (where r • r^ and n • nk) at an 
angle of from Snell's law for spherical stratification, one has: 


(lbw 


nfc r^ cos 8 ii • constant 


(187) 


If the angle Yjk due to refractive bending is computed for only 
a very small layer of atmosphere between the heights r * r^ and 


r - r k , the following assumptions are Justified: 

( 1 ) » constant, 

( 2 ) the index of refraction n is close to unity, and 

( 3 ) Ah - r^ - rj « rj (Ah is an infinitesimal layer 
of height), 

Using these assumptions along with Equations (186) and (187), Weiabrod 
and Andersonderive an expression for Yjk 1° terms of the angle 
6 and the refractivity N (see equation (168) for its definition).? 
for the details the reader is referred to their paper. Their expres¬ 
sion for Yjk through a small layer of atmosphere is: 

Uk ‘ 500 (t^.; + M» 6 k ) < " lllir * di *n ,) 

The total bending through the atmosphere is the sum of the 


( 188 ) 
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individual contributions. Thus, 


“ H. . - N. 

i-1 i „ 

y = Zu CAn TZ - a -TT-(milliradians) 

i=1 500 (tan p^ + tan p ) 


(189) 


From Equation (187), the value of cos p is given as: 

i 


cos p^ = 


"l-l Vl «■ ?!■! 

n. r. 

1 x 


Vi cos h-i 

n. 

l 


1 - 


Ah 


R + h. 
e l 


( 190 ) 


and it follows: 


J 1 - cos^ p^ 

tan P i = cos p. ” (191) 

1 

Equations (189), (190), and (191) determine the value of y 
where p is the value of the angle of incidence of the ray at a 
height tr = N(tr) is the refractivity at height tr obtained 
from the atmospheric model, and Ah is the increment of height used 
in Equations (189) and (190). Note that h^ = h + A* and that 
h^ is the computed height of the vehicle obtainable from the program. 

At the station h = h Q , which is the height of the station, and p^ 
is set equal to the computed elevation, E^. 

Once y has been found, consideration of the geometry in Figure 
27 leads to the evaluation of 6, the ’’error" in the elevation angle. 

Thus, we note from Figure 27: 


e = y - ( a - P) (192) 

Using Snell's law for spherical stratification and the application of 
the law of sines to triangle S0Q results in: 


cos P 


"o 


nr 


cos p 


0 


r Q = f3 0 = r cos a 

Combining Equations (193) and (194) yields: 

Q _ HO 

cos * n cos a = cos (<* ~ (a - j3)) 


(193) 

(194) 


(195) 
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Figure 27. Geometry of Ray Path Used to Obtain <5 


106 








Since a - (3 is a small angle, a small angle approximation and an 
approximation for is used^ 15 ^ in Equation (195) to obtain 

an expression for a - (3 in terms of refractivity at the station 
(Nq) and at the vehicle (N): 

a - 6 * (N 0 - K) 10* 6 cot 6 096) 

Using Equations (192) and (196), the value of c is determined, 

Applying the law of sines to triangle SOP results in: 


r 0 

Combining Equations 
tan 


cos (Bo - S) - r coi ((Jt + t) - O 

(194) and (197) gives an expression 

6 . tan a + .0 - cos g)_ 

•in c + cos c tan a - tan 8 q 


for 


(197) 


( 197 ) 


Finally, using snail angle approximations for 6 and 


£ 

tan 

a + 

e 2 1 

2 

* " e + 

tan 

a - 

tan f'o 


( 199 ) 


This value of 5 can now be used in Equation (176) to find the 
vector p 1 , which determines all angle related observations includin 
the effects of refraction. Due to the approximations used in the 
analysis, below elevations of about 5° the errors in the propagation 
corrections amount to a few percent. Because of this fact, propaga¬ 
tion corrections due to refraction should be computed only when the 
elevation of the vehicle exceeds 5 # . 


Refraction Effects on Range and Range Pxate 

Because of signal retardation caused by the refractive gradiirt 
of the atmosphere, the round-trip range computation of Equation (158) 
includes a correction term, Ap. The signal retardation caused by r»n 
infinitesimal layer of thickness dr is given by: 
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d T m - -i*1 esc 6 dr 
v V c J 


- ft -1) 


esc B 


dr 


- N x IO -6 & dr 
c 

The effect on range between layers bounded by heights r - rj and 
r - r^ is given by: 


(21 .<') 


r k 

r 


^pjk 


r k 

r 


cdT 


N x 10~ 6 CSC 6 dr 

* 

r J r J 

Using methods similar to the computation of retractive bending, 
Weisbrod and Anderson*" arrive at the expression for the refractive 

round-trip range correction, 2Ap, due to the passage of the ray 
through t\ie entire atmosphere: 


( 2 < 1 ) 


2Ac - 2 x io" 6 


y | Nl_l + I ( l 'i ” f-1 ^ 

sin Bi-i + sir. 


i-1 


+ IO -6 


1 + Tf 


l^i-i ” I 

^ sin Bi-i + sin 8^ 

i-irrH 1 1 


C' 


t. 


where 


is the refractivity in the i^ layer, 
h^ is the height of the i** 1 layer, 

6^ is obtained from equation (190), 
fj is the transmitted (or up) frequency, 
f 2 is the received (or down) frequency. 


The indices in the first summation refer to layers in the troposphere, 
and those in the second refer to layers in the ionosphere. The values 
of - N(h^) are obtained from the models of refractive index. 

Due to refractive bending, the range-rate measurement Of Equa¬ 
tion (159) also includes a correction term, Ap. This correction 
term arises because the program initially computes the component of 
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vehicle velocity along the vector p, i.e., along the straight line 
SP of Figure 27. The measured range-rate is the component of vehi¬ 
cle velocity along the tangent to the ray path at the vehicle, i.e., 
along the straight line TP of Figure 27. 

The velocity of the vehicle and the station is available from the 
program. If one denotes the component of vehicle velocity along the 
line SP by v x , and the component of vehicle velocity perpendicular 
to SP by Vy, then the range rate computation without the correc¬ 
tion for refraction (i.e., Ap - 0) is: 



The measured component of velocity along the tangent to the ray path 
at the vehicle as seen from Figure 27 is 


p m * v x cos (y - 6) + v y sin (y - 6) 
The correction to range-rate, Ap, is given by: 

AP - P m * Pc 


A factor should also be included to account for the round-trip fre¬ 
quency dependent effects in the ionosphere. After making a small angle 
approximation the final equation used for Ap, the one-way correc¬ 


tion to range-rate, is: 


Vy (Y - <5) 


where 



vy is the component of vehicle-station velocity perpendicular 
to the vector p, and lying in the plane which includes 
the center of the Earth, 

f\ is the transmitted frequency, 

f 2 is the received frequency. 


(203) 

(204) 

(205) 

(206) 
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SECTION III. 


USER'S GUIDE 

INTRODUCTION 

The information in this section is intended as a guide for those 
desiring to use SPACE-A. SPACE-A generates a vehicle ephemeris and 
associated observations from ground stations and has four modes of 
operation. They are: 

(1) Computation of vehicle ephemeris only; 

(2) Generation of observations from ground stations; 

(3) Same as (2) with the writing of the observations onto 
tape in a format suitable for processing by SPACE-B; 

(4) Same as (2) with the following included as standard 
outputs; 

(a) the acquisition time, 

(b) total time in sight, and 

(c) the option of the printing time of polar and 
meridian crossings. 

DESCRIPTION OF INPUT DATA 

All input data is read at the beginning of each run. The format 
and arrangement of the data is designed to make the data prepara¬ 
tion as convenient as possible. It is seldom necessary to read all 
the data; in fact, most runs require a minimum of cards. Many vari¬ 
ables have standard values which are preset by the program and should 
be satisfactory for most cases. 

The input data is divided into three parts. The first consists 
of a single card, called the basic input card, which must be included 
in each run. The second and third parts are designated group I and 
group II inputs, respectively. Each group is divided into several 
sections, each of which must be preceded by a card containing the 
section number, in integer format, in columns four and five. Any 
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section may be omitted from the input deck. Parts of sections, how¬ 
ever, must not be omitted. See Figure 28 for the setup of a sample 
input deck. 

The basic input card, which precedes each case, contains four 
variables. The first pertains to the standard values which may be 
set in place of the group II data and will be explained later. The 
second, MDE, indicates the function to be performed (i.e., trajectory 
computation, observables, etc.). The third indicates the desired pre¬ 
cision and determines the set of standard values used for integration. 
The last specified Encke or Cowell integration. 

The group I inputs consists of several sections. Any of these 
sections may be omitted from the input deck if the variables contained 
in the section are not needed for the run or in the case of stacked 
cases, the values from the last case are to be used again. The des¬ 
cription of each variable given in the data summary (see Table IV) is 
sufficient to understand most of the inputs. Therefore, only a brief 
discussion of the variables given in Table IV is presented here. 

In the group I input data, if the variable KLM2 of Section 2 is 
0 or 1, the initial conditions are assumed to be referred to the base 
date system or true system of date, respectively. 

Sections 7 and 9 deal with observations from observing systems 
on vehicles and with powered flight and have been omitted since they 
are not presently used. The variables in Section 8 control the out¬ 
put and are discussed in DESCRIPTION OF OUTPUT INFORMATION. 

The end of the group I input data is indicated by a card con¬ 
taining a 10 in columns four and five. This card must always be In¬ 

cluded in the data deck. After reading this card the program advances 
to the group II inputs. 

Before the program reads the group II input data, the val»e of 
the first variable contained on the basic input card, KSTDRD, is 

tested. If this variable is negative, standard values of the vari- 
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Figure 28. Sample Input Deck 






















































ables in the group II input data are computed from stored information 
and the program than begins to read the group II inputs to augment or 
replace any of the standard values. Table V gives a list of the stan¬ 
dard values used by the program. If KSTDRD is not negative, none of 
the standard values is set and all values should be read in unless 
values from a preceding stacked case are used. The end of group II 
inputs is indicated by reading a card with a 20 in integer format in 
columns four and five. 
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Table IV. 

Summary of Input Data 
Basic Input Card 


COLUMN NAME FORMAT 


DESCRIPTION 


2-5 KSTDRD 15 

6-10 MDE 15 


11-15 PRECIS F5.0 


16-20 CEPID F5.0 


< 0 Will want at least some 
standard inputs from 
group II 

^ 0 No standard. All values 

will be read in. 

- 1 Trajectory computation. 

* 2 Observable computations 
without simulated data. 

■ 3 Observable computations 

with simulated data. 

■ 4 Prediction mode 

■ 1 Low precision level. 

■ 2 Intermediate precision 

level. 

» 3 High precision level 

■ 0 f l Encke integration 

■ - 1 Cowell integration 
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Group I Input! 


SECTION CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

1 1 

2-72 

ITITLE(1-12) 

12A6 

Title 

2 

2-5 

NYEARP 

15 

Year of initial 
conditions 


6-10 

DAYS 

F5.0 

Day of initial 
conditions. 


11-15 

HRS 

F5.0 

Hour of initial 
day. 


16-20 

HMIN 

F5.0 

Minute of initial 
hour. 


21740 

SEC 

E20.16 

Seconds of initial 
minute. 

3 

1-24 

TMAX 

E24.16 

Time of run in 
hours. 

2 1 

2-5 

KLM 

15 

Indicator for units 
of position and 
velocity vector. 

- 1 KM,KM/SEC 

- 2 ER,ER/HR 

- 3 FT,FT/SEC 

- 4 MI,MI/HR 

- 5 NM,NM/HR 

- 6 NM,FT/SEC 

- 7 AU, AU /HR 


6-10 

KLM1 

15 

Indicator for coor¬ 
dinate system of 


input vectors. 

- 1 Cartesian coord. 

WE - 0 

■ 2 Cartesian coord. 
Compute WE 

- 3 Geodetic long, 

lat, alt, 
V e .V n ,V h ; WE - 0 

- 4 Geodetic long, 

lat, alt, 

Ve.V n ,V h ; 

Compute WE 
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SECTION CARD 


Groujs^J^JLnjjutis 

(continued) 

COLUMNS NAME FORMAT 


11-15 KLM2 15 


16-20 KSNAP 15 


21-25 KLM3 15 


DESCRIPTION 

■ 5 Geodetic long, 

lat, alt, 

|V|, Y, az; 

WE - 0 

- 6 Geodetic long, 

lat, alt, 

IVI, Y» az; 
compute WE 

* 7 Geocentric ra, 

decl, alt, 

v rat v d»Vh» 

WE - 0 

- 8 Geocentric ra, 

decl, alt, 
Vra,V d ,V h ; 
compute WE 

■ 9 Geocentric ra, 

decl, alt, 

|v| , 6, az; 

WE - 0 

* 10 geocentric ra, 

decl, alt, 

|v|, 6, az; 
compute WE 

Indicator for nuta¬ 
tion and precession 
of input vector. 

- 0 no 

- 1 yes 

Indicator for nuta¬ 
tion and precession 
of vectors during 
run. 

■ 0 no 

* 1 yes 

Indicator for libra- 
tion of input vectors. 

- 0 no 

- 1 yes 
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Group I Inputs 

(continued) 


SECTION 

CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 



26-30 

KLIBR 

15 

Indicator for libra- 
tion of vectors 
during run 

■ 0 no 

■ 1 yes 



31-35 

MR REF 

15 

Indicator for ini¬ 
tial reference body 






- 1 Earth 

- 2 Sun 

- 3 Moon 

- 4 Mars 

- 5 Venus 

- 6 Jupiter 

- 7 Saturn 


2 

1-24 

RCIN(l) 

E24.16 

Initial position 
vector 



25-48 

(2) 

E24.16 

See KLM and KLM1 
for units and type. 



49-72 

(3) 

E24.16 



3 

1-24 

RCIN(l) 

E24.16 

Initial velocity 
vector. 



25-48 

(2) 

E24.16 

See KLM and KLM1 
for units and type. 



49-72 

(3) 

E24.16 


3 

1 

2-5 

PASFX 

F5.0 

Total number of 
passes. 

4 

1 

2-5 

KS2BY 

15 

Indicator for two- 
body integration 

- 0 no 

- 1 yes 



6-10 

KSPLT 

15 

Indicator for in¬ 
clusion of planetary 
perturbations 

- 0 no 
» 1 yes 
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SECTION 


Group I Inputs 
(continued) 

CARD COLUMNS NAME FORMAT 

11-15 KSOBL 15 


16-20 KSDRG 15 


21-25 KSRAP 15 


26-30 KSDRGM 15 


31-35 KSDRGV 15 


36-40 KSMNOB 15 


41-45 KRF 15 


DESCRIPTION 

Indicator for inclu¬ 
sion of oblateness 
perturbation# 

- 0 no 

■ 1 yes 

Indicator for inclu¬ 
sion of Earth drag 
perturbation. 

■ 0 no 

■ 1 yes 

Indicator for inclu¬ 
sion of radiation 
pressure perturba¬ 
tion. 

■ 0 no 

- 1 yes 

Indicator for inclu¬ 
sion of Mars drag 
perturbation. 

- 0 no 

- 1 yes 

Indicator for inclu¬ 
sion of Venus drag 
perturbation. 

■ 0 no 

- 1 yes 

Indicator for inclu¬ 
sion of Moon oblate¬ 
ness perturbation. 

■ 0 no 

- 1 yes 

Indicator for inclu¬ 
sion of refraction 
effects. 

■ 0 no 

- 1 yes 
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Group I Inputs 


(continued) 

SECTION CARD COLUMNS NAME FORMAT DESCRIPTION 

46-50 KECLPS 15 Indicator for inclu¬ 

sion of eclipse in¬ 
formation print. 


- 0 no 

- 1 yes 


5 1 

2-5 

MAXSTA 

15 

Total number of 
stations used in 

run. 

2 

2-5 

K 

15 

Station number, 


10-15 

STANM(L) 

A6 

Station name. 


20-30 

TYPE(L) 

4X111 

A(1)+l,E2*A(2) 


+l.E4*A(3)+i.E6*A(4) 
+1,E8*K where 
A ■ observation 
types used by 
station K in 
ascending order, 
L - packed station 
number 


3 

1-24 

STALT(K) 

E24.16 

Latitude of station 

K, 


25-48 

STALN (K) 

E24.16 

Longitude of station 

K. 


49-72 

STAHT(K) 

E24.16 

Altitude of station 
K. 

4 . 

1-36 

RRATE(1,L) 

E36.16' > 

C Repetition rates 

1 (hrs). 


37-72 

RRATE(2,L) 

E36.16 

\ \ For each observation 

5 

1-36 

RRATE(3,L) 

E36.16 

vType 


37-72 

RRATE(4,L) 

E36.16> 


6 

1-18 

TDELAY(1,L) 

E18.8 

C Times in hrs before 


19-36 

TDELAY(2,L) 

E18.8 

J which each observa- 


37-54 

TDELAY(3,L) 

E18.8 

S tion is not to be 


55-72 

TDELAY(4,L) 

E18.8 

computed. 
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Group I Inputs 

(continued) 


SECTION CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

7 

1-18 

FUP(K) 

E18.8 

Station transmit 
freq. (MHz). 


19-36 

FDOWN(K) 

E18.8 

Station receive 
freq. (MIIz). 

8 

1-24 

STA0R(NCDST+1) 

E24.16 

AEE station 
rotation angle. 


25-48 

STA0R(NCDST+2) 

E24.16 

AEV station 
rotation angle. 


49-72 

STAOR(NCDST+3) 

E24.16 

AEN station 
rotation angle. 

9 

1-24 

STAOR(NCDST+4) 

E24.16 

AU 

^station lo¬ 
cation 


25-48 

STAOR(NCDST+5) 

E24.16 

-< 

errors 
caused by 


49-72 

STAOR(NCDST+6) 

E24.16 

AW 

geodetic 
jiet error. 

10 

1-24 

STAOR(NCDST+7) 

E24,16 

NO; refractiv- 
vity at sea 
level. 


25-48 

STAOR(NCDST+8) 

E24.16 

H; troposphere 
scale fact. 


49-72 

STAOR(NCDST+9) 

E24.16 

P0; max. elec¬ 
tron density. 

11 

1-24 

STAOR(NCDST+10) 

E24.16 

HO; 

of 

loser limit 
ionosphere. 


25-48 

STAOR(NCDST+11) 

E24.16 

HM; Ht of PO 
(KM) 


49-72 

STAOR(NCDST+12) 

E24.16 

- open - 

12 

1-24 

STAOR(NCDST+13) 

E24.16 

AT 

for timing. 


25-48 

STAOR(NCDST+14) 

E24.16 

Bias added for 
obser.A(l) 


49-72 

STAOR(NCDST+15) 

E24.16 

Bias added for 
obser.A(2). 
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Group I Inputs 
(continued) 


SECTION CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

13 

1-24 

STAOR(NCDST+16) 

E24.16 

Bias added for 
obser*A(3) 


25-48 

STAOR(NCDST+17) 

E24.16 

Bias added for 
obser•A(4) 


REPEAT 

CARDS 2-13 FOR EACH 

STATION 


6 1 

1-18 

DAREAS 

E18.8 

Effective sur¬ 
face area of 
vehicle per¬ 
taining to 
drag (ft 2 ). 


19-36 

PAREAS 

E18.8 

Effective sur¬ 
face area per¬ 
taining to 
radiation pres¬ 
sure (ft 2 ). 


37-54 

SPADD(6) 

E18.8 

Mass of vehicle 
(lb-masses) 

7 

OPEN 




8 1 

2-5 

I UNIT 

15 

Indicator for 
output units 
(see KLM) 


6-55 

IPSEC(I),I“1.10 

1015 

Indicator for 
suppression of 
each of 10 print 
sections. 

2 

1-18 

19-36 

DTP I 

DTSUP 

E18.8 

E18.8 

Print portion 
(hrs) and sup¬ 
press portion 
(hrs), of total 
print period. 


37-54 

PRATE 

E18.8 

Interval with¬ 
in DTPI at which 
to print. 


9 OPEN 

10 END OF GROUP I INPUTS 
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Group II Inputs 


SECTION CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

1 1-7 

1-72 

(DT3(I,J), 
1-1,3) J-1,7 

3E24.12 

Integration in¬ 
tervals for each 
of 7 working 
bodies for near, 
medium and far 
reference. 

2 1 

1-72 

R1(1-3) 

3E24.12 

R1 and R2 are 

2 

1-72 

R1(4-6) 

3E24.12 

distances in 

E.R. for each of 

3 

1-72 

R1(7),R2(1-2) 

3E24.12 

7 working bodies 

4 

1-72 

R2(3-5) 

3E24.12 

for switching 
from near to 

5 

1-48 

R2(6-7) 

3E24.12 

medium and me¬ 
dium to far in¬ 
tegration inter¬ 
vals* 

3 1 

1-24 

RT1 

E24.12 

Values used as 


25-48 

RT2 

E24.12 

tolerances in 



rectification 
criteria. 



4 1 

1-18 

DH1 

E18.8 

Troposphere 
integration 
step size (KM)* 


19-36 

DH2 

E18.8 

Ionosphere in¬ 
tegration step 
size (KM). 


37-54 

H2 

E18.8 

Upper limit 
troposphere 
(KM). 


55-72 

H4 

E18.8 

Upper limit 
ionosphere (KM)• 

5 1 

2-5 

KOBLAT 

15 

Number of ob¬ 
lateness coeffi¬ 
cient terms* 

2 

2-5 

N 

15 

N index. 


6-10 

M 

15 

M index. 

3 

1-24 

SORC1 

E24.12 

Value of C 
coefficient. 
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Group II Inputs 



(continued) 



SECTION CARD COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

25-48 

S0RC2 

E24.12 

Value of S 
coefficient 


REPEAT CARDS 2-3 UNTIL KOBLAT VALUES HAVE BEEN READ IN. 


6 

1 

2-5 

MBMAX 

15 

Number of working 
bodies. 



6-10 

KWBMU(I) 

1*1,MBMAX 

15 

Indices of 
working bodies. 


2 

1-72 

TPMAT8(1) 
1-1,MBMAX 

3E24.12 

Gravitational 
constants of 
working bodies. 


REPEAT 

CARD 2 

FOR EACH VALUE 

OVER 3N NEEDED. 

7 

1 

1-24 

DYN(48) 

E24.12 

Solar flux in 
10 -22 w/m 2 -Hz 
at 10.7 cm. 



25-48 

DYN(49) 

E24.12 

Open 

8 

1 

1-72 

DYN(51-53) 

3E24.12 

Coefficients 
for iunar ob¬ 
lateness 
(kg • km 2 ). 

9 

1 

1-24 

COMB(l) 

E24.12 

Velocity of 
light. 

10 

1 

1-72 

PRNT3(l-3) 

3E24.12 

Print intervals 
(hrs) for near 
medium and far 
reference. 

11 

1 

1-24 

EMIN 

E24.12 

Minimum value 
of elevation 
angle (radians). 

12 

1 

1-18 

RTO 

E18.8 

Ratio of 

Nordsieck inte¬ 
gration interval 
to that in Runge- 
Kutta. 
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SECTION 

13 


14 


15 


16 

17 

18 

19 

20 


Group II Inputs 
(continued) 


CARD 

COLUMNS 

NAME 

1 

1-24 

DSPL 

1 

1-72 

RATEC(1-3,1) 

2 

1-72 

RATEC(l-3,2) 

1-10 

1-72 

XMACH(I) 

1-1,40 

11-20 

1-72 

CDT(I) 

1-1,40 


FORMAT 

DESCRIPTION 

E24.12 

Special inte¬ 
gration inter¬ 
val in A4 mode 
to obtain ac¬ 
quisition time 
(hrs). 

3E24.12 

Rotation vector 
used in Mars 
drag computa¬ 
tions. 

3E24.12 

Rotation vector 
used in Venus 
drag computa¬ 
tions. 

4E18.8 

Mach number 
table. 

4E18.8 

Drag coefficient 
table 


OPEN 

OPEN 

OPEN 

OPEN 

END OF GROUP II INPUTS 
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Table V. 


Standard Values for Group II Inputs 


Section 1. 

DT3 ARRAY (3 * 7 

DT3(1,I); 1-1,3,4,5,6,7 
DT3(2,I) ; 1-1,3,4,5,6,7 
DT3(3,I); 1-1,3,4,5,6,7 
DT3(1,2) - 4. 

DT3(2,2) - 6. 

DT3(3,2) - 10. 


array of Integration Intervals) 
.1953125 x io- 2 
.15625 x 10" 1 
.125 


Note: The above values are always used for Encke integration and for 
Cowell when low precision is specified. When Cowell and inter¬ 
mediate precision is specified, each value in the DT3 array is 
divided by 2; when Cowell and high precision is specified each 
value is divided by 4, The precision desired is controlled by 
the variable CEPID on the Basic Input Card, 


Section 2, 

Rl ARRAY (1x7 array of reference body change 

criteria) 

Rl (I) - 4 • RADII(I) 

R2 (I) - 10 • RADII(I) 

where RADII(I) is the radius of the I*h working body (Earth, Sun, 
Moon, Venus, Mars, Jupiter, Saturn), 

Section 3, 

RTl - 1. x 10" 6 
RT2 - 1. x 10" 6 
If low precision is specified 
RTl - 1 x 10" 4 
RT2 - 1 x IQ" 4 
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Table V. 

(continued) 


Section 4. 

DH1 - 1. 

DH2 - 10. 

H2 - 40. 

H4 - 2,000. 

Section 5. 


KOBLAT - 3 


C 20 " “ 

1082.3 * 

10~ 6 

s 20 

c 30 " 

X 

CO 

• 

CM 

10" 6 

s 30 

c 4 0 * 

1.8 « 

10 -6 

S 40 


Section 6. 

KBMAX - 7 

KWBMU(I) - I; I - 1,2,3,4,5,6,7 


TPMAT8(1) 

- 

19.909416 

(2) 

- 

6629965.8 

(3) 

- 

.24478289 

(4) 

- 

16.1983009 

(5) 

- 

2.14364682 

(6) 

- 

6338.16258 

(7) 

- 

1897.36734 


Section 7. 

DYN(48) - 200. 

Section 8. 

DYN(51) - .88746 * 10 29 kg • km 2 

DYN(52) - .88764 « 10 29 kg • kn 2 

DYN(53) - .88807 * 10 29 kg * ka 2 


128 


Table V. 
(continued) 


Section 9. 

C0MB(1) - 169210.58004 

Section 10. 

PRNT3(1) - .25 
PRNT3(2) - .5 
PRNT3(3) - 1. 

Section 11. 

EMIN - (5. degrees; expressed in radians) 

Section 12. 

RTO - 3 


Section 13. 

DSPL - 1. x 10" 3 


Section 14. 

RATEV(1,1) 
RATEV(2,1) 
RATEV(3,1) 
RATEV(1,2) 
RATEV(2,2) 
RATEV(3,2) 


0.0 


0.0 

0.255175469 

0.0 

0.0 


0.0 
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Table V. 
(concluded) 


Section 15. 

MACH TABLE ARRAY (XMACH(I), I - 1,40) 


.5 

1.5 

25. 

96. 

0.0 

1.85 

30. 

97. 

.5 

2.0 

50. 

98. 

.75 

2.7 

70. 

99. 

.85 

3.4 

90. 

100. 

.9 

4.2 

91. 

101. 

1.04 

5.6 

92. 

102. 

1.1 

6.75 

93. 

103. 

1.2 

8.5 

94. 

104. 

1.3 

15. 

95. 

105. 

DRAG COEFFICIENT TABLE (CDT(I) 

.8 

1.4 

1.14 

1.14 

.8 

1.38 

1.14 

1.14 

.82 

1.36 

1.14 

1.14 

.92 

1.285 

1.14 

1.14 

.99 

1.23 

1.14 

1.14 

1.04 

1.19 

1.14 

1.14 

1.18 

1.16 

1.14 

1.14 

1.21 

1.155 ■ 

1.14 

1.14 

1.26 

1.15 

1.14 

1.14 

1.3 

1.45 

1.14 

1.14 


40) 
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DESCRIPTION OF OUTPUT INFORMATION 

The input quantities which determine the form of the output of 
the program are given in Section 8 of the Group I inputs in Table IV. 
The choice of output units are determined by IUNIT. The choices of 
units with corresponding values of IUNIT are given below. 

Value of IUNIT OUTPUT UNITS 

1 KM,KM/SEC 

2 ER,ER/HR 

3 FT,FT/SEC 

4 MI,MI/HR 

5 NM,NM/HR 

6 MM,FT/SEC 

7 AU,AU/HR 

Time may be divided into periods of print and suppression of 
print by setting the input variables DTPI and DTSUP. A period of 
printing (DTPI) Is executed followed by a period of suppression (DTSUP) 
and the cycle is then repeated. 

The variable PRATE, is used to determine the print interval with¬ 
in DTPI. Since conditions occur where the print routine is entered 
more than once with the same time, printing is normally only executed 
upon the first entry* However, by using a negative PRATE, printing 
will occur each time the print subroutine is entered during the print 
interval determined by DTPI. 

The user has the option of choosing the information he would like 
printed by setting IPSEC(I) * 1, where I is the output section number 
given in Table VI, that contains certain information. By reading 
a zero for IPSEC(I), the printing of the I** 1 section is surpressed. 
There are ten sections of output information as shown in Table VI. See 
Figure 29 for a sample listing with various sections labeled. 
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Table VI. 

Ten Sections of Output Information 


Section 1. 

Program time in hours since launch; 

Current program date in hours, minutes, seconds, 
day and year; 

Planetary body which is the current reference. 

Section 2. 

Station number; 

Observation types which station handles. 

Section 3. 

Components and magnitudes of position and velocity 
vectors in units specified by the user. 

Section 4. 

Components and magnitudes of perturbations in 
position, velocity and acceleration in units 
specified by the user. 

Section 5. 

Components and magnitudes of vectors between 
vehicle and each planetary body in specified 
units. 

Section 6. 

Components and magnitudes of vectors between 
the Earth and each of the other planetary bodies 
in specified units. 

Section 7. 

Components and magnitudes of vectors between 
the Sun and each of the planetary bodies 
in specified units. 
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Table VI. 
(concluded) 


Section 8, 

Right ascension; 

Declination; 

Greenwich hour angle; 

Longitude and latitude of Earth subsatellite point 
Geocentric altitude, azimuth and elevation; 
Geodetic azimuth and elevation. 

Section 9. 

True, eccentric and mean anomalies; 

Semi-major axis; 

Eccentricity; 

Inclination; 

Argument of perigee; 

Perigee passage time; 

Right ascension of ascending node; 

Period; 

Mean motion; 

Perigee and apogee heights; 

Unit vector to perigee; 

Unit angular momentum vector. 

Section 10. 

Station name; 

Program time in hours; 

Station location; 

Observation values. 
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PROGRAM DECK SETUP 

This section describes the deck setup for running the A-mode of 
SPACE on the 7030 computer. The program is executed under the control 
of MCP (Master Control Program). The binary cards containing the com¬ 
piled routines are kept on a magnetic tape at the MITRE 7030 computer 
facility and only a few control cards and the data deck are needed by 
the user. 

Those cards shown in Figure 30 that have a B or T punched in 
column 1 are control cards for the MCP. They are described in greater 
detail in the MITRE 7030 Facility Manual, 

JOB Card 

The job card is always the first card of an input deck. This 
card contains necessary information for accounting and no job can be 
run without it. The fields on the JOB card are variable in length 
and separated by commas. 


Format 
1 10 

B JOB, 1 JOBID*,’NAME*,’PROJECT*,'DEPT*,’MAXTIME 1 ,'DEST* 


Field Name 

Description 

B 

Punched in column 1 

JOB, 

Punched in columns 10-13 

'JOBID' 

1 to 7 character Job identi¬ 
fication 

'NAME' 

Name of programmer (up to 13 
characters) 

'PROJECT’ 

Project number (3 to 9 
characters) 

'DEPT' 

Department number (3 
characters) 

'MAXTIME' 

Maximum running time in 
minutee (up to 3 characters) 

'DEST' 

Destination of output 
(up to 6 characters) 


136 





Example: 

1 10 

B JOB,SPACE,BOND J, 007, D84 , 13, A141 

TYPE Card 

This card defines the type of job (GO), and the proper monitor 
system (FORTRAN) to the MCP to execute this program. 


Format: 

1 10 

B TYPE,GO,FORTRAN 

ADDIQ Card 

This card contains the library number of the program tape that 
contains the binary decks of the SPACE program. Since modifications 
to the program may be made in the future, the user should periodically 
check the current tape number with Department D-84. 


Format: 

10 

ADDIO,PLBXXXX 


1 

BLIB1 


Field Name 
BLIB1 

ADD10,PL BXXXX 


SUBTYPE,BIN Card 

This card tells the MCP that 


Description 

Punched in columns 1-3 

Beginning in column 10, 
the XXXX indicate the 
field where the tape 
number is to be punched. 


e cards that follow are binary. 


Format: 

1 10 

T SUBTYPE,BIN 

TMAIN Card 

This card results in the construction of a dummy main program 
which calls the executive routine of SPACE (EXECA). EXECA is compiled 
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as a subroutine and the binary cards are contained on the library tape* 
This card is the equivalent to the Fortran routine* 

CALL EXECA 

RETURN 
END 


Format: 

1 10 

T MAIN,EXECA 

Processed FIOD Deck 

Since the programmer cannot choose the input-output channels to 
be used (they are assigned by MCP at execution time), the FIOD deck 
specifies the I/O channel requirements of the program* The processed 
FIOD deck is the result of compiling the following subprogram; 


1 

10 


T 

SUBTYPE,FIOD 


B2 

IOD,$READER 


B3 

IOD,SPRINTER 


B9 

IOD,TAPE,,,,,, 

.SAVE 

B8 

IOD. TAPE.. 

.SAVE 

B 

REEL,PLBDL680 



END 



Logical units 2 and 3 are the system input and output units res¬ 
pectively. Logical unit 9 is used for the output of the observation 
tape generated when that option is exercised* Logical unit 8 is used 
for the input of the ephemeris tape containing the positions of the 
Moon and planets* The reel card specifies that the ephemeris tape, 
whose label is PLBL680, is to be mounted* 

SUBTYPE.DATA Card 

This card tells the MCP that the following cards contain the 
data for the program* 
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Formats 


1 10 

T SUBTYPE,BIN 

Input Data Deck 

This is the data deck described in the Description of Input Data. 



Figure 30. Deck Setup 
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SECTION IV. 
PROGRAMMING 


GENERAL FLOW OF PROGRAM 

The general flow of the program is shown in Figure 31* The user 
does not need to specify the number of cases to be run, as the program 
returns to read the data for the next case until all data cards are 
exhausted. 

The program structure and linkage of the program between basic 
routines is shown in Figure 32. It is to be noted, however, that the 
figure does not indicate every linkage but only the major ones. For 
instance, in some cases it is necessary for OBSERA (which usually only 
computes observations) to call the integration package. 

The routines listed under the heading of General Routines, in 
Figure 32 are used by many of the major routines. The linkage of 
these routines is not included as this would defeat the purpose of 
the figure, which is to indicate the position of the major routines 
in the overall structure of the program. 

DESCRIPTION OF SUBROUTINES 

The following part of the document gives a brief description of 
the subroutines used by SPACE-A. Each subroutine description gives 
the purpose and a brief description of the subroutine, as well as the 
other subroutines which call and are called by the given subroutine. 
For understanding the details of each subroutine one ought to refer 
to the pertinent portions of Section II and Section III dealing with 
the function of that subroutine as well as referring to the computer 
listing. 

It should be noted that, as of the writing of this document, 
certain subroutines have not been tested* These non-operational sub¬ 
routines are: LUNOBL, MVDRAG, OBD, PFINIT, PFLGHT, STACUL. These 
subroutines deal with lunar oblateness gravitational perturbations. 
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BECIN 


RETURN TO 
FIND NEXT 
ACTIVITY 
TIME 



RETURN 
TO READ 
DATA FOR 
NEXT CASE 


Figure 31. General Flaw of Program 
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*872 


EXECA 

(EXECUTIVE ROUTINE) 


INPUT A 

(READS DATA) 


MAIN A 

(CONTROLS MAIN 
FLOW OF PROGRAM) 


XFQRM 

(READS 8 CONVERTS 
INITIAL CONDITIONS) 


TIMNGA 

(FIND NEXT 
ACTIVITY TIMS 


OBSER A 

(COMPUTES 
OBSERVATIONS) 


PFINIT 

(INITIALIZE 
POWER ED FLIGHT! 


STACUL 

(CHECK FOR 
OCCULATION) 


PRINT A 

(PRINT SELECTED 
INFORMATION) 


COWELL 



INTEGRATION 
ROUTINES 



OBD (ON-BOARD OBSERVATIONS) 


* ATIM (MODE 4 COMPUTATION) 


STAPOS (STATION VECTOR) 


ENCKE 


♦ MODEL A (INDEX OF REFRACTION) 


CITGRA (CONTROLS INTEGRATION) 
CCHREF ( CHANGE REFERENCE BODY) 
CRSTRE ( SAVE/RESTORE VECTORS) 
CINT (INTEGRATION) 


b| t )' 


GENERAL ROUTINES 

DDOT (COMPUTE A B > ) 

DMTML (COMPUTE Mfe|/IA 

DOMUD (CHECK FOR ERROR) 

EPHEM (READ EPHEMERIS TAP£L[ 

FIX (UNPACK WORD) 

NUTPRE (COMPUTE NUTAT./ 
PREC.) 

SERVCE (COMPUTE AXB, 

ICI) 


EITGRA (CONTROLS INTEGRATION) 
ECHREF ( CHANGE REFERENCE BODY) 
ERSTRE (SAVE/RESTORE VECTORS) 
EINT (INTEGRATION) 

KEPLER (NOMINAL/2-BODY ORBIT) 
RECT (RECTIFICATION OF ORBIT) 


I 


PER IV 

(COMPUTES ACCELERATIONS) 




I 


OBLDRG 

(OBLATENESS 
AND DRAG) 


LUNOBL 

(LUNAR 
OBLATENESS) 


MVDRAG 

(MARS/VENUS! 


DENSTY 

(COMPUTES 
DENSITY 
OF AIR) 


ECLIPSE 


PFLGHT 

(CHECK a PRINT 


(POWERED 

ECL. INFO.) 


FLIGHT) 


Figure 32. Program Structure 
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Mars and Venus drag model, on-board observation * $ perturbations due 
to vehicle thrust, and observations dealing with occultation. 

EXECA 

Purpose : 

This is the executive routine for SPACE-A* 

Description : 

The program alternates between calling INPUTA and MAINA until 
all cases have been exhausted. 

Subroutines Called : INPUTA, MAINA 

ATIM 

Purpose: 

This subroutine is used in the visibility computation mode (mode 
4). It determines the time a vehicle comes within sight of a given 
ground station or the time of polar baseline crossing (the time that 
the vehicle crosses the east-west vertical plane through the zenith) 
and meridian crossing* 

Description : 

When the vehicle is in sight, the l and m direction cosines 
are tested and used to interpolate for the times of polar and meri¬ 
dian crossings* When the acquisition time is to be determined, the 
integration routines CINT and EINT or the two-body routine KEPLER is 
employed to predict ahead until the vehicle comes into view. 

Called By : OBSERA 

Subroutines Called : CINT, CRSTRE, EINT, ERSTRE, KEPLER 

CCHREF 

Purpose: 

This subroutine tests criteria for changing the reference body 
when the Cowell integrator is used* Even though a change of reference 
body in the Cowell method is not necessary, it is utilised in the 
program in the same manner as in the Encke integrator* 
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Description ! 

Criteria for changing the reference body based on regions of 
influence are tested. 


Called By : CITGRA 

Subroutines Called : DDOT, EPHEM, SERVCE 


CINT (EINT) 

Purpose ; 

This subroutine perforins the Integration when the Cowell formu¬ 
lation is employed. 

Description ! 

The Gill modification of Runge-Kutta is used for short-term inte¬ 
gration and as a starting procedure for the Nordsieck long-term inte¬ 
gration. 

Called By : ATIM, CITGRA, MAINA, OBSERA 
Subroutines Called ! CCHREF, CINT, CRSTRE 

CITGRA 

Purpose : 

This subroutine serves as the sub-main program governing calls 
to the Integration routines in the Cowell method. 

Description : 

If not in powered flight, the program checks for a change of re¬ 
ference bodies. The delta of integration is determined by the dis¬ 
tance of the vehicle from the refetence body and integration is per¬ 
formed up to the next activity time by calling the integration routine. 
Called By : MAINA 

^ubroutinj^^g^ecH CCHREF, CINT, CRSTRE 

CRSTRE (ICR) 

Purpose : 

This subroutine saves or restbres time, position, and velocity 
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and acceleration of the vehicle at any designated time for the Cowell 
integrator. 

Description ; 

Depending on the value of the argument (IRC) the information is 
saved (IRC ^ 1) or restored (IRC > 1), 

Called By ; ATIM, CITGRA, MAINA, OBSERA 

DDOT(A.B) 

Purpose: 

This function computes the dot product of two vectors. 
Description ; 

The program uses the formula 

DDOT - A(l)* B(1) + A(2)* B(2) + A(3)* B(3) 

Called By ; CCHREF, DERIV, ECHREF, LUNOBL, MVDRAG, OBD f 
OBLDRG, OBSERA, PRINTA, RECT, STACUL 

DENSTY(XI, X2, X3, DENP t DENEQ) 

Purpose ; 

This subroutine evaluates the density of air in the upper 
atmosphere. 

Description ; 

Tables of the logarithm of p • C av obtained from the Harris- 
Priester data are stored in this routine. Interpolation is performed 
to calculate log(p • C aV ) (DENEQ) at the equator as a function of the 
altitude of the vehicle (XI), solar flux (X2), and solar time (X3)• 
Interpolation is also performed to calculate log(p • Cav)(DENP) at 
the poles, as a function of the altitude and solar flux. 

Called By : OBLDRG 

DERIV 

Purpose: 

This subroutine evaluates the acceleration terms for the Encke 
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and Cowell integrators. 

Description : 

This subroutine computes the planetary perturbations and the solar 
radiation pressure perturbations as presented in the theory section of 
this manual. Earth oblateness perturbations and drag are computed by 
calling OBLDRG. Provisions are included for computing powered flight 


accelerations. 

Called By : CINT, EINT 

Subroutines Called ; DDOT, DOMUD, ECLIPS, EPHEM, KEPLER, 

LUNOBL, MVDRAG, OBLDRG, PFLGHT, SERVCE 

DMTML (A, B, C t I t J, K, L, M, N, IAC, JAB, KBC, IFLAG) 

Purpose: 


This subroutine multiplies two matrices of any size (up to 26 * 
26), or two matrices in which the second is transposed. 

Description : 

The matrices A, B, and C are dimensioned I * J, K * L, and 
L * M respectively. 


When IFLAG “0, 1, multiplication is done by rows of A so 
that A can be overwritten if desired (i.e., the argument C in the 
calling sequence may actually be the same as A). 

When IFLAG “2, 3, multiplication is done by columns of B and 
the result is stored in B. 

The following table summarizes the results of the program for 
various values of IFLAG. 


IFLAG 

0 

1 

2 

3 


COMPUTE 

T 

A • B 
A • B 

T 

A . B 
A • B 


STORE RESULT IN 

C 

C 

B 

B 


IAC is the number of rows of A and C to be used, 

JAB is the number of columns of A and rows of B 
(or B t ) to be used, 
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IBC is the number of columns of B (or B?) and 
C to be used. 

Called By: LUNOBL, NUTPRE, OBLDRG, OBSERA, 

STACUL, STAPOS, XFORM 

DOMUD (TEST) 

Purpose : 

This subroutine checks for errors. 

Description : 

The Overflow and Divide Check indicators are checked. If either 

is on, and the argument TEST is not equal to zero, TEST is assumed to 

be a BCD word and the message "ERROR IN 'TEST'" is printed. 

Called By: DERIV, KEPLER, NUTPRE, OBLDRG, OBSERA, 

PRINTA, RECT, STAPOS, XFORM 

ECHREF 

Purpose : 

This subroutine tests criteria for changing the reference body 
when the Encke integrator is used. 

Description : 

Criteria for changing the reference body based on regions of 
influence are tested. 

Called By : EITGRA 

Subroutines Called : DDOT, EPHEM, KEPLER, SERVCE 

ECLIPS (J, K) 

Purpose : 

This subroutine checks for a change in the illumination of the 
vehicle and prints out an appropriate message when the vehicle 
changes between the three states: umbra, penumbra, or sunlight. 
Description : 

This first argument in the calling sequence indicates whether 
the reference body is the Earth (J»l), the Moon (J*2), or some other 
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body (J>2), The second argument indicates whether the vehicle is in 
the umbra (K-l), penumbra (K-2), or sunlight (K*3)» Values of the 
second argument are saved from one entry to the next and hence a 
simple comparison is made for determining a change in vehicle status. 

If a change has occurred! a message is printed containing the refer¬ 
ence body, the time of change and the entered state (umbra, penumbra, 
sunlight). 

Called By : DERIV 

EINT (IENT) 

Purpose ; 

This subroutine performs the integration when the Encke method 
is employed. 

Description : 

The Gill modification of Runge-Kutta is used for short-term in¬ 
tegration and as a starting procedure for the Nordsieck long-term 
integration. 

Called By : DERIV 

Subroutines Called : ATIM, EITGRA, OBSERA 

EITGRA 

Purpose : 

This subroutine serves as the sub-main program governing calls 
to the integration routines in the Encke method. 

Description : 

If not in powered flight, the program checks for a change of re¬ 
ference bodies. The delta of integration is determined by the distance 
of the vehicle from the reference body and integration is performed up 
to the next activity time by calling the integration routine. Fre¬ 
quent checks are also made on the magnitude of the deviations from the 
nominal orbit. If the magnitudes increase beyond specified limits, a 
rectification of the orbit is accomplished by calling RECT. 
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Called By: MAINA 


Subroutines Called : ECHREF, EINT, ERSTRE, KEPLER, RECT 

EPHEM 

Purpose: 

This subroutine evaluates the position and velocity vector* of 
each of the seven bodies with respect to the reference body* 
Description ! 

Tabular planetary positions are read from an ephemcris tape and 
Everett’s interpolation formula is fitted to six points* It is eval¬ 
uated to find the position vector and its derivative is evaluated to 
find the velocity. 

Called By ; CCKREF, DERIV, ECHREF, STaCUL 
Subroutines Called : SERVCE 

ERSTRE (IRC) 

Purpose : 

This subroutine saves or restores time, position, velocity, and 
acceleration of the vehicle at any designated time for the Encke in¬ 
tegrator. 

Description : 

Depending on the value of the argument (IRC) the information is 
saved (IRC ^ 1) or restored (IRC > 1). 

Called By : ATIM, EITGRA, MAINA, OBSERA 

FIX (KTEMP, ITEMP, KNAME) 

Purpose : 

This subroutine unpacks a word into five separate words. 
Description : 

The argument KTEMP is assumed to contain a number of the form 
VVWWXXYYZZ. The word is unpacked and stored as follows: 

KNAME - W 
ITEMP(1) - ZZ 
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ITEMP(2) - YY 
ITEMP(3) - XX 
ITEMP(4) - WW 
Called By : OBSERA 

INPUTA 

Purpose : 

This subroutine reeds all data necessary for one run. 

Description : 

The subroutine Initializes necessary data and reads In sections 
desired. 

Called By : EXECA 
Subroutines Called : XFORM 

KEPLER 

Purpose : 

This subroutine computes the two-body position and velocity 
vectors. 

Description : 

A Newton Raphson scheme Is used to determine the differential ec¬ 
centric anomaly. After convergence, the tworbody position end velo¬ 
city vectors are evalueted. The subroutine uses Herrick's method. 
Called By : ATIM, DERIV, ECHREF, EITGRA, OBSERA, STACUL 
Subroutines Called : SERVCE, DOMUD 

LUNOBL(K) 

Purpose : 

This subroutine computes the acceleration due to lunar oblateness. 
Optionally, It can compute the libratlon and effect of the eerth's 
J 20 term. 

Description: 

When K ■ 1, the libratlon matrix is computed and then precessed 
and nutated. 
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Whan K - 2 , tha Earth’s J 20 oblatanaaa tarn is calculatad. 
When K - 3, both lunar and Zarth oblatanaaa (J20 tars) ara 
computed. 

Called By : DER1V, OBD 

Subroutines Called : DDOT, DWTML, NUTPRE, SERVCE 

MAINA 

Purpose : 

This subroutine handles the main flow of the program. 

Descrj^Hti^ii: 

The routine calls subroutines which determine the next activity 
time, govern the integration, compute the observations, and print the 
chosen information. 

Called By : EXECA 

Subroutines Called : CITGRA, CRSTR1, ZITGRA, ERSTRE, OBSIRA, 

PFINIT, PRINTA, *ECT, SERVCE, STACUL, 

TIMNGA 


MODELA(K) 

Purpose : 

This subroutine computes the index of refraction for the tropo¬ 
sphere or ionosphere. 

Description : 

The index of refraction is computed with the tropospheric model 
when K ■ 1 and the ionospheric model when K ■ 2. 

Called By: OBSERA 


MVP RAG 


Purpose : 

This subroutine computes the perturbations due to drag in tha 
atmospheres of Mars or Venue. 

Description : 

The coefficient of dreg is determined by interpolation from given 
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tables as a function of altitude and the velocity of the vehicle. 

Called By : DERIV 

Subroutines Called : DDOT, SERVCE 

NUTPRE(K) 

Purpose : 

This subroutine computes the expressions for determining the 
Greenwich hour angle, as well as the rotational matrices between co¬ 
ordinate systems. 

Description : 

When K ■ 1, the expressions used In the nutation matrix and the 
llbratlon matrix are computed. 

When K ■ 2, the rotation matrix through the Greenwich hour 
angle Is computed. 

When K ■ 3, the precession-nutation matrix Is computed. 

Called By : LUNOBL, OBLDRG, PRINTA, STAPOS, XFORM 
Subroutines Celled I DMTML, DOMUD 

OBD 

Purpose : 

This subroutine computes observations from on-board Instrumen¬ 
tation. 

Description : 

The subroutine uses present position of vehicle with respect to 
reference bodies, landmarks, or ground stations to determine obser¬ 
vations. 

Called By : OBSERA 

Subroutines Called : DDOT, LUNOBL, SERVCE, STAPOS 

OBLDRG 

Purpose : 

This subroutine computes the oblateness and air drag perturbations 
due to the Earth and Its atmosphere. 
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Description; 


Oblateness is computed in a flexible algorithm which can be used 
for zonal, tesseral, or sectorial harmonics of any order. The drag 
perturbation is computed from the Harris-Priester tables (at high alti¬ 
tudes) and the U. S. Standard Atmosphere (at low altitude*). 

Called By ; DERIV 

Subroutines Called ; DDOT, DENSTY, DMTML, DO MUD, NUTPRE, SERVCE 

OBSERA 

Purpose: 

This subroutine computes the observables as seen from a given 
ground station. 

Description : 

Given present vehicle location and ground station location, both 
referenced to inertial space, the observed values are computed. When 
specified, refraction corrections are included. When in mode 4, the 
total time in sight at the station is also computed. 

Called By ; MAINA 

Subroutines Called : ATIM, CINT, CRSTRE, DDOT, DMTML, DOMUD, 

EINT, ERSTRE, FIX, KEPLER, MODELA, OBD, 

SERVCE, STAPOS 

PFINIT 

Purpose : 

This subroutine performs initialization procedures at the entry 
of a burn period. 

Description : 

The subroutine determines coefficients of six Chebyshev polyno¬ 
mials which are valid for the first burn period. 

Called By; MAINA 
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PFLCHT 

Purpose: 


This routine computes the effect of a thrust perturbation in the 
Gncke method. 

Description : 

The thrust acceleration vector at the initial burn time, the vehi¬ 
cle mass and mass rate are converted to a set of coefficients for 6 
Chebyshev polynomials. These coefficients are determined in PFINIT, 

The subroutine, PFLCHT, uses these coefficients to describe the powered 
flight trajectory as a function of time. Effectively, the subroutine 
replaces KEPLER in Gncke's method. Integration of the equations of 
motion continues exactly as if powered flight was not Involved. 

Called By : DERIV 
Subroutines Called; servce 

PRINTA 

Purpose: 

This subroutine prints out current trajectory information and 
observations. 

Description : 

Depending on whether the time entered corresponds to a print ‘time, 
the program checks each of the ten elements in the IPSEC array. If the 
value is > 0, the corresponding section is printed. 

To determine whether it is a print time, the subroutine first 
checks to see whether the present time is within the print portion 

(DTPI) of a total print period (TAU). If not, no printing is done. 

If so, it next checks the value of the print interval within DTPI 
(PRATE). If it is negative, it always prints. If it is positive, 
and it is the first time into the present print period, it prints, 

otherwise no printing is done. When MDE ■ 1 or 4, printing occurs 

at each entry to the routine. 
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Called By i MAINA 

Subroutines Called : DDOT, DOMUD, NUTPRE, SERVCE 


RECT 

Purpose ; 

This subroutine computes the parameters pertinent to a rectifi¬ 
cation in Encke’s method. 

Description : 

The two-body position and velocity vectors are equated to the 
true position and veolcity vectors. In addition, these components are 
saved. Perturbations in position and velocity are set equal to zero, 
and elements used by the KEPLER subroutine are computed. 

Called By : EITGRA, MAINA 

Subroutines Called : DDOT, DOMUD 

SERVCE (A, B, C, I) 

Purpose : 

This subroutine computes the cross product of two vectors as well 
as the magnitude, magnitude squared and magnitude cubed of a vector. 

Description : 

When 1*1, the cross product of A and B is computed and 
the result stored in C(l), C(2), and C(3), and then continues as 
when 1*2. 

When 1-2, the magnitude cubed, magnitude, and magnitude 
squared is stored in C(4), C(5), and C(6), respectively. 

Called By: CCHREF, DERIV, ECHREF, EPHEM, KEPLER, LUNOBL, 

MAINA, MVDRAG, OBD, OBLDRG, OBSERA, PFLGHT, 

PRINTA, STACUL 


STACUL 

Purpose : 

This subroutine determines the time of occultation (YCOM(23)). 
One of two types of occultation may be considered, vehicle or star 
occultation* 
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Description : 

A Newton-Raphson iteration is used on the two-body solution to 
determine the time of occultation. When considering star occultation, 
the occulting body is the reference body, except in the Earth-Moon 
system, in which either the Moon or the Earth can do the occulting. 
Vehicle occultation occurs only in Moon reference. Up to three dif¬ 
ferent ground stations may be considered for occultation. 

Called By : MAINA 

Subroutines Called : DDOT, DMTML, EPHEM, KEPLER, SERVCE, STAPOS 

STAPOS 

Purpose : 

This subroutine computes the station position and velocity vectors. 

Description : 

Given the geodetic latitude, longitude, and altitude of the 
station, the coordinates of the station in the Greenwich coordinate 
system are computed. A rotation through the Greenwich hour angle then 
brings the coordinates into the inertial system. The velocity of the 
station is computed from the rotational speed of the earth. 

Called By : OBD, OBSERA, STACUL 

Subroutines Called : DMTML, DOMUD, NUTPRE 

TIMNGA 

Purpose : 

This subroutine determines the next activity time (i.e,, the next 
time at which an observation is to be computed or a printout of the 
trajectory is to occur). 

Description : 

Logic is set up for establishing an array of times of interest 
from which the earliest time is selected. Flags are set when the 
time selected is the final time or when e time is repeated. 
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Called By : MAINA 


XFORH 

Purpose : 

This routine accepts the input information concerning the vehi¬ 
cle’s initial conditions, and transforms them into the proper units 
and coordinate system. 

Description : 

After the initial conditions and associated flags are read in 
(Section 2 of the group I inputs), the program converts them into 
the position and velocity vectors in the units of the Earth radii, 
and Earth radii/hour. Then, if desired, the vectors are transformed 
into the base date system. 

Called By : INPUTA 

Subroutines Called : DHTML, DOMUD, NUTPRE 
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